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Josephson  junction  (JJ)  integrated  circuits  (ICs)  are  capable  of  operating  at  clock 
frequencies  from  1  to  100  GHz.  At  these  frequencies,  analysis  of  signal  propagation 
delay,  crosstalk,  dispersions,  radiation,  and  reflections  must  be  included  to  determine 
proper  response  of  the  circuit.  Much  effort  is  required  in  simulating  high  frequency 
behavior,  where  the  cross-sectional  dimensions  of  conductors  are  comparable  to  the 
signal  wavelength,  with  conventional  circuit  simulation  methods  as  SPICE.  A  simulation 
method  capable  of  modeling  high-frequency  behavior  by  solving  Maxwell’s  curl 
equations,  the  finite-difference  transmission  line  matrix  method  (FD-TLM),  is  modified 
to  model  JJ  logic  circuits  and  provide  simultaneous  time-domain  three-dimensional  full- 
wave  electromagnetic  field  and  JJ  device  analysis.  The  FD-TLM  method  is  further 
extended  to  model  superconducting  quantum  interference  devices  (SQUIDs). 

Techniques  for  simulation  and  simulation  results  for  a  Josephson  Atto-Weber  switch 
(JAWS),  a  two-junction  superconducting  quantum  interference  device  (SQUID),  and  a 
modified  variable  threshold  logic  (MVTL)  gate  are  provided.  Interconnection  lengths  are 
kept  intentionally  short  so  the  FD-TLM  simulations  can  be  validated  by  conventional, 
low  frequency,  quasi-static  analysis.  The  general  behavior  observed  in  FD-TLM 
simulation  and  good  agreement  with  quasi-static  conventional  circuit  simulation  validate 


the  FD-TLM  method. 


FD-TLM  SIMULATION  OF  JOSEPHSON  JUNCTION  LOGIC 

CIRCUITS 


Christopher  G.  Sentelle,  MS 
University  of  Nebraska,  1995 


Advisor:  Robert  H.  Voelker 


Josephson  junction  (JJ)  integrated  circuits  (ICs)  are  capable  of  operating 
at  clock  frequencies  from  1  to  100  GHz.  At  these  frequencies,  analysis  of  signal 
propagation  delay,  crosstalk,  dispersion,  radiation,  and  reflections  must  be  included 
to  determine  proper  response  of  the  circuit.  Much  effort  is  required  in  simulating 
high  frequency  behavior,  where  the  cross-sectional  dimensions  of  conductors  are  com¬ 
parable  to  the  signal  wavelength,  with  conventional  circuit  simulation  methods  as 
SPICE.  A  simulation  method  capable  of  modeling  high-frequency  behavior  by  solv¬ 
ing  Maxwell’s  curl  equations,  the  finite-difference  transmission  line  matrix  method 
(FD-TLM),  is  modified  to  model  JJ  logic  circuits  and  provide  simultaneous  time- 
domain  three-dimensional  full- wave  electromagnetic  field  and  J  J  device  analysis.  The 
FD-TLM  method  is  further  extended  to  model  superconducting  quantum  interference 
devices  (SQUIDs).  Techniques  for  simulation  and  simulation  results  for  a  Josephson 
Atto-Weber  switch  (JAWS),  a  two-junction  superconducting  quantum  interference 
device  (SQUID),  and  a  modified  variable  threshold  logic  (MVTL)  gate  are  provided. 
Interconnection  lengths  are  kept  intentionally  short  so  the  FD-TLM  simulations  can 
be  validated  by  conventional,  low  frequency,  quasi-static  analysis.  The  general  be¬ 
havior  observed  in  FD-TLM  simulation  and  good  agreement  with  quasi-static  con¬ 
ventional  circuit  simulation  validate  the  FD-TLM  method. 


References 

[1]  R.  H.  Voelker  and  R.  J.  Lomax,  “Three-dimensional  simulation  of  integrated 
circuits  using  a  variable  mesh  TLM  method,”  Microwave  and  Opt.  Technol.  Lett., 
vol.  2,  pp.  125-127,  Apr.  1989. 

[2]  M.  Kamon,  C.  Smithhisler,  and  J.  White,  “FastHenry  User’s  Guide,”  Research 
Laboratory  of  Electronics,  Dept,  of  Elect.  Engr.  and  Comp.  Sci,  Massachusetts 
Institute  of  Technology,  Cambridge,  MA,  Sept.  13,  1994. 

[3]  B.  Johnson,  et  al.,  “SPICE3  version  3el  User’s  Manual,”  Dept,  of  Elect.  Engr. 
and  Comp.  Sci.,  Univ.  of  California,  Berkeley,  CA,  Apr.  1,  1991. 

[4]  B.  W.  Char,  et  al.,  First  Leaves:  A  Tutorial  Introduction  to  Maple  V.  New  York: 
Springer- Verlag,  1992. 

[5]  W.  H.  Press,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P.  Flannery,  Numerical 
Recipes  in  C ,  2nd  ed.  New  York:  Cambridge  University  Press,  1992. 

[6]  S.  Hasuo  and  T.  Imamura,  “Digital  Logic  Circuits,”  Proceedings  of  the  IEEE., 
vol.  77,  no.  8,  pp.  1177-1190,  Aug.  1989. 

[7]  S.  Takada,  “Progress  on  Josephson  Computer,”  Nonlinear  Superconductive  Elec¬ 
tronics  and  Josephson  Devices.,  Plenum  Press:  New  York,  1991. 


121 


122 


[8]  H.  Jackel,  W.  Bachtold,  “Computer-Simulation  for  Digital  Josephson  Devices 
and  Circuits,”  Circuit  Analysis,  Simulation  and  Design.  North-Holi  and:  Elsevier 
Science  Publishers,  1986 

[9]  J.  G.  Rollins,  “Numerical  simulator  for  superconducting  integrated  circuits,” 
IEEE  Trans.  Computer-Aided  Design,  vol.  10,  pp.  245-251,  Feb.  1991. 

[10]  J.  Clarke,  “Principles  and  Applications  of  SQUIDs,”  Proceedings  of  the  IEEE., 
vol.  77,  no.  8,  pp.  1208-1221,  Aug.  1989. 

[11]  T.  P.  Orlando  and  K.  A.  Delin,  Foundations  of  Applied  Superconductivity. 
Addison- Wesley  Publishing  Company:  New  York,  1991. 

[12]  R.  H.  Voelker  and  C.  G.  Sentelle,  “FD-TLM  Modeling  of  Josephson  Junction 
Circuits,”  IEEE  Microwave  and  Guided  Wave  Lett.,  vol.  5,  pp.  137-138,  May 
1995. 

[13]  T.  VanDuzer  and  C.  W.  Turner,  Principles  of  Superconductive  Devices  and  Cir¬ 
cuits.  New  York:  Elsevier,  1981. 

[14]  S.  T.  Ruggiero  and  D.  A.  Rudman,  Superconducting  Devices.  Boston:  Academic 
Press,  1990. 

[15]  _ ,“A  finite-difference  transmission  line  matrix  method  incorporating  a  nonlinear 

device  model,”  IEEE  Trans.  Microwave  Theory  Tech.,  vol.  38,  pp.  302-312,  Mar. 
1990. 

[16]  R.  H.  Voelker  and  C.  G.  Sentelle,  “Graphical  modeling  interface  for  FD-TLM 
electromagnetic  field  simulation  of  high-speed  pulse  propagation  in  active  struc¬ 
tures,”  in  Dig.  1994  IEEE  AP-S  Inti.  Symp.,  Seattle,  WA,  June  1994,  pp.  1120- 
1123. 


123 

[17]  G.  Costabile,  S.  Pagano,  N.  F.  Pedersen,  and  M.  Russo,  Nonlinear  Superconduc¬ 
tive  Electronics  and  Josephson  Devices.  New  York:  Plenum  Press,  1991. 

[18]  M.  A.  Megahed,  “Nonlinear  Analysis  of  Microwave  Superconductor  Devices  Us¬ 
ing  Full-Wave  Electromagnetic  Model,”  IEEE  Transactions  on  Microwave  The¬ 
ory  and  Techniques .,  vol.  43,  no.  11,  pp.  2590-2598,  Nov.  1995. 

[19]  C.  G.  Sentelle  and  R.  H.  Voelker,  “FD-TLM  Electromagnetic  Field  Simulation 
of  High-Speed  Josephson  Junction  Digital  Logic  Gates,”  IEEE  Transactions  on 
Microwave  Theory  and  Techniques.,  accepted  for  publication,  Dec.,  1995. 


FD-TLM  SIMULATION  OF  JOSEPHSON  JUNCTION  LOGIC 

CIRCUITS 


by 

Christopher  G.  Sentelle 


A  THESIS 


Presented  to  the  Faculty  of 
The  Graduate  College  at  the  University  of  Nebraska 
In  Partial  Fulfillment  of  Requirements 
For  the  Degree  of  Master  of  Science 


Major:  Electrical  Engineering 


Under  the  Supervision  of  Professor  Robert  H.  Voelker 


Lincoln,  Nebraska 


December,  1995 


FD-TLM  SIMULATION  OF  JOSEPHSON  JUNCTION  LOGIC 

CIRCUITS 

Christopher  G.  Sentelle,  MS 
University  of  Nebraska,  1995 


Advisor:  Robert  H.  Voelker 


Josephson  junction  (JJ)  integrated  circuits  (ICs)  are  capable  of  operating 
at  clock  frequencies  from  1  to  100  GHz.  At  these  frequencies,  analysis  of  signal 
propagation  delay,  crosstalk,  dispersion,  radiation,  and  reflections  must  be  included 
to  determine  proper  response  of  the  circuit.  Much  effort  is  required  in  simulating 
high  frequency  behavior,  where  the  cross-sectional  dimensions  of  conductors  are  com¬ 
parable  to  the  signal  wavelength,  with  conventional  circuit  simulation  methods  as 
SPICE.  A  simulation  method  capable  of  modeling  high-frequency  behavior  by  solv¬ 
ing  Maxwell’s  curl  equations,  the  finite-difference  transmission  line  matrix  method 
(FD-TLM),  is  modified  to  model  JJ  logic  circuits  and  provide  simultaneous  time- 
domain  three-dimensional  full-wave  electromagnetic  field  and  JJ  device  analysis.  The 
FD-TLM  method  is  further  extended  to  model  superconducting  quantum  interference 
devices  (SQUIDs).  Techniques  for  simulation  and  simulation  results  for  a  Josephson 
Atto-Weber  switch  (JAWS),  a  two-junction  superconducting  quantum  interference 
device  (SQUID),  and  a  modified  variable  threshold  logic  (MVTL)  gate  are  provided. 
Interconnection  lengths  are  kept  intentionally  short  so  the  FD-TLM  simulations  can 
be  validated  by  conventional,  low  frequency,  quasi-static  analysis.  The  general  be¬ 
havior  observed  in  FD-TLM  simulation  and  good  agreement  with  quasi-static  con¬ 
ventional  circuit  simulation  validate  the  FD-TLM  method. 


ACKNOWLEDGEMENTS 


I  would  like  to  give  special  thanks  to  Dr.  Rober  H.  Voelker  for  his  support  and 
encouragement  during  the  term  of  this  project.  I  would  also  like  to  thank  Dr.  Voelker 
for  his  understanding  during  my  difficult  semester  in  the  fall  of  1994.  I  would  like  to 
give  thanks  to  Dr.  Alexander  and  the  Center  for  Electro-Optics  for  financial  support 
during  this  project.  I  would  like  to  thank  my  parents  for  their  encouragement  and 
understanding  in  my  life.  Finally,  I  would  like  to  give  special  thanks  to  my  wife  for 
her  unending  understanding,  support,  and  encouragment. 


Table  of  Contents 

1  Introduction  1 

2  Background  and  Rationale  4 

3  The  FD-TLM  Method  for  Josephson  Junction  Simulation  8 

3.1  Modeling  Linear  Devices .  12 

3.2  Modeling  Nonlinear  Devices . 14 

4  Modeling  the  Josephson  Junction  16 

4.1  Characteristics  of  the  Josephson  Junction .  16 

4.2  Implementation  of  the  Josephson  Junction  in  the  FD-TLM  method  .  21 

5  Resistively  Coupled  Josephson  Logic  Circuits  31 

5.1  Introduction .  31 

5.2  The  JAWS  circuit .  32 

5.2.1  FD-TLM  implementation  .  34 

5.2.2  Validation  of  FD-TLM  results .  38 

6  Magnetically  Coupled  Josephson  Logic  Circuits  42 

6.1  Introduction .  42 

6.2  DC  SQUID  Simulation .  42 


i 


ii 

6.2.1  FD-TLM  implementation  .  44 

6.2.2  Validation  of  FD-TLM  results .  45 

6.3  MVTL  Simulation  .  47 

6.3.1  FD-TLM  implementation  .  48 

6.3.2  Validation  of  the  FD-TLM  results .  50 

7  Conclusion  53 

A  FORTRAN  Code  for  JJ  Implementation  56 

B  JAWS  Simulation  Data  65 

B.l  FD-TLM  Data  Set .  65 

B.2  FastHenry  Data  Set .  69 

B.3  MAPLE  V  Results .  70 

B. 4  C  Code  .  73 

C  2  JJ  DC  SQUID  Simulation  Data  79 

C. l  FD-TLM  Data  Set .  79 

C.2  FastHenry  Data  Set . 86 

C.3  MAPLE  V  Results .  88 

C. 4  C  Code  .  .  .  . .  96 

D  MVTL  Simulation  Data  102 

D. l  FD-TLM  Data  Set .  102 

D.2  FastHenry  Data  Set .  109 

D.3  MAPLE  V  Results .  Ill 

D.4  C  Code 


116 


r 


List  of  Figures 

3.1  Depiction  of  divided  three-dimensional  space .  10 

3.2  Each  cell  in  three-dimensional  space .  10 

3.3  The  E  field  node  .  11 

3.4  The  H  field  node .  11 

3.5  Equivalent  circuit  of  the  JJ  showing  current  source .  15 

4.1  DC  current- voltage  characteristics  of  the  JJ  .  17 

4.2  Equivalent  circuit  for  the  Josephson  junction .  18 

4.3  FD-TLM  simulation  for  a  single  Josephson  junction .  25 

4.4  Depiction  of  the  arrangement  of  JJs  within  the  SQUID  loop .  26 

4.5  Total  critical  current  versus  flux  for  the  SQUID .  28 

5.1  The  JAWS  logic  gate .  32 

5.2  Cross-section  of  layers  used  in  fabrication  of  JJ  ICs .  35 

5.3  Placement  of  the  Ey  node  in  the  tunneling  barrier .  35 

5.4  JAWS  layout  used  for  FD-TLM  simulation .  36 

5.5  Results  of  simulation  of  the  JAWS  logic  gate .  37 

5.6  MAPLE  is  used  to  solve  first  order  derivatives .  39 

6.1  Equivalent  circuit  for  the  2  JJ  DC  SQUID .  43 

6.2  Physical  layout  of  the  DC  SQUID .  44 

iii 


IV 

6.3  Results  of  simulation  of  the  DC  SQUID  logic  gate .  45 

6.4  Equivalent  Circuit  for  the  MVTL  logic  gate .  47 

6.5  Physical  layout  of  the  MVTL  circuit  used  in  FD-TLM  simulation.  .  .  49 

6.6  Results  of  simulation  of  the  MVTL  logic  gate .  50 


Chapter  1 
Introduction 


The  Josephson  junction  is  a  superconducting  tunneling  device  that  allows  fabrication 
of  logic  gates  with  much  lower  power  dissipation  and  much  faster  switching  speeds 
than  logic  gates  utilizing  silicon  or  GaAs  technology.  Because  Josephson  junctions 
are  capable  of  operating  at  much  higher  frequencies,  circuit  simulation  is  significantly 
more  complex.  Quasi-static  conventional  circuit  becomes  inaccurate  at  frequencies 
where  the  wavelength  of  the  signal  becomes  comparable  to  the  cross-sectional  dimen¬ 
sions  of  the  circuit  layout.  Thus,  a  full-wave  simulation  method  must  be  found  that 
is  capable  of  modeling  JJ  circuits  at  high  frequencies. 

The  full- wave  finite-difference  transmission  line  matrix  (FD-TLM)  method  [1], 
completely  models  the  time  evolution  of  electromagnetic  field  interaction  with  dif¬ 
ferent  media  and  devices  by  solving  Maxwell’s  curl  equations.  In  this  method,  three 
dimensional  space  is  divided  into  a  set  of  nodes  which  are  selectively  modified  to  model 
material  properties  such  as  permittivity,  permeability,  and  conductivity.  Nonlinear 
devices  are  modeled  in  the  method  by  altering  the  conductivity  and  permittivity  of 
the  nodes  representing  the  device  as  a  function  of  voltage  and  time. 

The  FD-TLM  method,  by  solving  Maxwell’s  curl  equations,  naturally  simulates 
the  operation  of  a  circuit  while  simultaneously  determining  high  frequency  effects. 
The  JJ  is  implemented  as  a  nonlinear  device  in  the  method  with  conductivity  and 
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permittivity  that  change  with  applied  voltage  and  time.  Substrates  and  interconnects 
are  modeled  by  specifying  the  values  of  permeability,  conductivity,  and  permittivity 
at  the  nodes  representing  them.  The  FD-TLM  method  can  model  JJ  logic  circuits  as 
well  as  the  substrate  layers  and  insulating  layers  used  to  fabricate  the  logic  circuit. 
Additionally,  there  is  no  need  to  extract  parasitic  component  values  since  these  are 
automatically  simulated  while  solving  Maxwell’s  equations.  The  FD-TLM  method  is 
discussed  in  Chapter  3. 

Many  logic  circuits  utilizing  JJs  use  a  property  of  JJs  that  arises  when  several 
JJs  are  placed  in  a  superconducting  loop.  That  is,  two  or  more  JJs  placed  in  the 
same  superconducting  loop  become  quantum  mechanically  coupled  with  a  relation¬ 
ship  controlled  by  the  presence  of  a  magnetic  flux  within  the  superconducting  loop. 
This  configuration  is  referred  to  as  a  Superconducting  Quantum  Interference  Device 

(SQUID). 

The  Josephson  junction  and  SQUID  are  implemented  in  the  FD-TLM  method  as 
described  in  Chapter  4.  Chapter  5  discusses  the  resistively  coupled  Josephson  logic 
gate  which  utilizes  the  characteristics  of  individual  JJs  while  Chapter  6  discusses  the 
magnetically  coupled  Josephson  logic  gate  which  utilizes  magnetic  coupling  between 
several  JJs  placed  in  a  superconducting  loop. 

The  results  from  the  FD-TLM  method  can  be  validated  by  observing  the  general 
logic  operation  obtained  from  the  simulation.  In  addition,  further  validation  is  ob¬ 
tained  by  performing  a  conventional  circuit  simulation.  A  program  FastHenry  [2]  is 
used  to  extract  inductance  values  from  the  physical  layout  of  the  JJ  logic  gate.  The 
inductance  values  are  required  to  accurately  simulate  JJ  circuits  when  using  a  conven¬ 
tional  circuit  simulation  method.  In  contrast,  the  FD-TLM  method  does  not  require 
a  separate  inductance  extraction  since  it  solves  Maxwell’s  curl  equations  rigorously. 

Since  methods  such  as  SPICE  [3]  do  not  currently  implement  the  JJ  device,  con- 
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ventional  circuit  simulation  proceeds  with  nodal  analysis  of  the  equivalent  circuit  for 
the  logic  gate  to  obtain  a  set  of  simultaneous  differential  equations.  MAPLE  V  [4] 
is  used  to  symbolically  solve  these  differential  equations  for  the  standard  form  where 
the  first-order  derivatives  are  on  the  left-hand  side  of  the  equations.  Then,  the  new 
set  of  differential  equations  are  solved  numerically  using  a  fifth-order  Runge-Kutta  [5] 
method,  providing  the  response  of  the  circuit  with  time.  Validation  of  the  FD-TLM 
method  for  each  type  of  logic  gate  is  discussed  in  Chapters  5  and  6. 

Chapter  7  ends  with  a  summary  discussion  of  the  results  obtained  from  the  simu¬ 
lation  of  all  logic  gates  along  with  a  statement  on  the  ability  of  the  FD-TLM  method 

to  simulate  circuits  incorporating  JJ  devices.  Possible  future  work  in  this  area  will 
be  discussed  as  well. 


Chapter  2 

Background  and  Rationale 


Theoretically,  superconducting  digital  logic  gates  should  be  capable  of  achieving  high 
speeds  with  low  power  consumption.  However,  early  superconducting  digital  circuits 
dating  back  to  the  1950s  utilized  cryotrons  and  were  slow  by  semiconductor  standards. 
Cryotron  digital  circuits  relied  on  the  physicall  change  between  a  superconducting  and 
normal  state  to  switch  between  a  logic  “1”  and  “0”  [6].  Resultant  switching  times 
were  a  few  microseconds.  With  superconducting  digital  technology  unable  to  compete 
with  semiconductor  technology,  further  research  was  abandoned. 

Research  into  the  use  of  the  Josephson  junction  (JJ)  for  superconducting  digital 
logic  circuits  began  in  the  mid-1960s.  The  Josephson  junction  utilizes  tunneling  of 
Cooper  pairs,  superconducting  transport  particles,  across  a  thin  tunneling  barrier 
consisting  of  a  non-superconductive  material.  Early  JJs  utilized  lead  alloy  supercon¬ 
ductors  with  Si02  as  tunneling  barriers.  The  tunneling  barrier  is  designed  to  allow 
tunneling  of  Cooper  pairs  but  also  allow  tunneling  of  normal  electrons.  As  is  typical 
for  all  superconductive  materials,  the  JJ  can  only  support  up  to  a  maximum  super¬ 
conducting  current  density.  Once  this  current  density  is  reached,  the  device  switches 
to  conduction  of  normal  electrons  through  the  thin  barrier  rather  than  conduction 
of  superconducting  Cooper  pairs.  In  other  words,  the  device  switches  between  a 
superconducting  and  normal  conducting  state.  This  is  done  without,  for  example, 
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changing  the  temperature  of  the  device  meaning  that  much  higher  switching  speeds 
can  be  obtained  than  with  the  cryotron. 

Although  JJ  logic  circuits  exhibit  high  speed,  low  power  performance,  silicon 
ICs  are  still  used  in  most  applications  today  for  several  reasons.  Research  into  JJs 
and  superconducting  ICs  declined  after  a  short  burst  in  the  1960s  due  to  the  low, 
liquid  Helium,  temperatures  required  for  operation,  the  inability  to  precisely  control 
parameters  of  the  JJ  in  fabrication,  and  the  dramatic  changesin  parameters  of  the  JJ 
with  extended  storage  and  temperature.  As  a  consequence  of  the  impracticality  of 
use  of  JJs  in  IC  circuits  at  the  time,  IBM  discontinued  research  into  JJ  technology 
in  1983. 

With  the  advent  of  Type  II  superconductors  capable  of  operating  at  liquid  ni¬ 
trogen  temperatures,  JJs  recently  became  more  practical  for  digital  logic  ICs.  In 
addition  to  operating  at  much  higher  temperatures,  JJs  using  superconductive  ma¬ 
terials  such  as  YBaCuO  and  Nb  are  much  more  stable  than  their  predecessors  in 
that  their  parameters  such  as  critical  current  do  not  change  significantly  over  time. 
Furthermore,  it  is  easier  to  reproduce  JJs  with  similar  electrical  properties  on  an  IC 
with  the  newer  superconductive  materials. 

Unfortunately,  while  several  JJs  can  be  fabricated  with  the  same  properties,  it  is 
still  impossible  to  fabricate  ICs  with  more  than  a  few  thousand  JJs  all  having  similar 
properties.  Therefore,  while  small  scale  integration  is  feasible,  it  is  still  impossible  to 
create  superconducting  ICs  with  the  level  of  integration  used  in  silicon  technology. 
Furthermore,  even  though  Josephson  junction  logic  gates  consume  far  less  power 
and  have  much  higher  switching  speeds  than  their  counterparts  in  GaAs  or  silicon 
technology,  the  performance  of  GaAS  and  silicon  technology  is  rapidly  approaching 
that  of  JJ  technology  while  remaining  more  economical. 

Despite  the  stumbling  blocks,  research  is  still  ongoing  in  the  design  of  even  faster 
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JJ  logic  gates.  Many  JJs  being  fabricated  today  have  switching  times  as  low  as  2.5 
ps.  Furthermore,  several  4-bit  microprocessor  architectures  utilizing  JJs  have  been 
fabricated  and  tested.  Of  particular  note,  in  1988,  a  team  of  researchers  developed 
a  4-bit  microprocessor  using  2.5  pm  Nb  technology  [7].  The  processor  consisted  of  a 
64-bit  Non-Destructive  Read-Out  (NDRO)  RAM,  an  ALU,  and  control  circuits  using 
a  total  of  1841  Josephson  junction  gates.  The  microprocessor  operated  at  a  maximum 
clock  frequency  of  770  MHz. 

Research  is  still  ongoing  in  developing  many  types  of  logic  gates  utilizing  the 
JJ.  Research  time  and  cost  can  be  minimized  given  an  ability  to  simulate  JJ  logic 
circuits  before  fabrication.  Accurate  circuit  simulation  reduces  the  amount  of  trial- 
and-error  fabrication.  In  simulating  JJs  that  can  switch  in  only  2.5  ps,  propagation 
delays,  dispersion,  reflections,  and  parasitic  capacitance  and  inductance,  must  all  be 
considered. 

Many  methods  have  been  used  to  simulate  JJ  [8]  circuits  including  a  recent  method 
where  SPICE  is  extended  to  model  JJs  [9].  Rollins  [9]  provides  a  starting  point  for 
implementing  JJ  circuits  in  other  simulation  methods.  However,  in  Rollins’  approach, 
SPICE  was  used  to  simulate  a  JJ  circuit  with  nanosecond  switching  times  and  was  not 
developed  to  ensure  proper  modeling  for  higher  clock  frequencies.  Additionally,  the 
SPICE  method  presented  by  Rollins  appears  to  omit  several  important  parasitics  in¬ 
cluding  mutual  inductance  between  the  control  line  and  superconducting  loop,  capac¬ 
itance  between  the  control  line  and  superconducting  loop,  and  capacitance  between 
the  circuit  and  the  superconducting  ground  plane.  However,  stray  capacitances  and 
inductances  can  be  extremely  important  as  illustrated  in  [10].  Here  it  is  emphasized 
that  parasitic  capacitance  between  the  control  line  and  superconducting  loop  can 
create  resonances  in  a  SQUID  circuit. 

In  summary,  the  SPICE  method,  while  simulating  the  characteristics  of  the  JJ 


itself  rather  well,  is  inept  at  modeling  JJ  logic  circuits  where  the  wavelength  of  the 
signals  become  comparable  to  the  cross-sectional  dimensions  of  conductors  in  the 
circuit.  For  instance,  with  the  small  voltage  signal,  approximately  2  mV,  used  in 
JJ  logic  circuits,  electromagnetic  noise  may  determine  whether  the  circuit  functions 
correctly.  A  method,  then,  must  be  used  that  is  capable  of  solving  Maxwell’s  equa¬ 
tions  to  provide  the  full  electromagnetic  response  of  the  circuit  and  its  surrounding 
packaging. 

Several  methods  exist  which  solve  Maxwell’s  curl  equations  for  electromagnetic 
fields.  Most  notably,  the  finite-difference  time-domain  (FD-TD)  method  which  solves 
discretizes  three  dimensional  space  into  a  set  of  nodes  and  solves  for  Maxwell’s  curl 
equations  in  time  and  space  using  the  “leap-frog”  method  for  integrating  the  differ¬ 
ential  equations.  The  finite-difference  transmission  line  matrix  (FD-TLM)  method 
works  in  a  similar  manner,  however,  this  method  solves  for  voltages  and  currents 
rather  than  electric  and  magnetic  fields  by  using  an  analogy  between  transmission 
line  equations  and  Maxwell’s,  curl  equations. 


Chapter  3 

The  FD-TLM  Method  for 
Josephson  Junction  Simulation 


Several  methods  exist  for  simulating  the  Josephson  junction  device.  SPICE  [9],  for 
example,  can  be  used  to  obtain  the  behavior  of  a  Josephson  junction  logic  circuit. 
Conventional  circuit  analysis  methods  such  as  SPICE  are  excellent  for  simulation  of 
circuits  including  high  speed  circuits  given  that  the  parasitics  such  as  inductance  and 
capacitance  are  extracted  from  the  physical  chip  layout  or  are  supplied  as  part  of  the 
circuit  for  simulation.  However,  extraction  of  parasitic  inductance  and  capacitance 
from  a  physical  chip  layout  can  be  inaccurate  and  may  only  be  correct  for  certain 
frequencies.  Phenomena  such  as  the  skin  effect,  lumped  versus  distributed  capaci¬ 
tance  and  inductance,  dispersion,  electromagnetic  noise,  and  resonance  within  the  IC 
package  cannot  be  readily  modeled  in  a  SPICE  simulation. 

The  FD-TD  and  FD-TLM  methods,  by  solving  Maxwell’s  curl  equations,  com¬ 
pletely  model  the  time-evolution  of  electromagnetic  field  interaction  with  different 
media  and  devices.  As  a  result,  high  speed  integrated  circuits  can  be  simulated  as 
with  the  SPICE  method  while,  in  addition,  simultaneously  modeling  distributed  par¬ 
asitic  capacitance  and  inductance,  dispersion,  electromagnetic  noise  and  reflections, 
and  cavity  resonance  given  the  physical  layout  of  materials  and  devices  within  the 
circuit. 
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While  the  FD-TD  method  solves  Maxwell’s  curl  equations  directly  using  a  dif¬ 
ferential  equation  integration  technique,  the  FD-TLM  method  relies  on  an  analogy 
between  Maxwell’s  curl  equations  and  the  equations  describing  transmission  lines. 
As  a  result,  the  FD-TLM  method  actually  solves  for  voltages  and  currents  that  are 
propagated  through  space  along  interconnecting  transmission  lines  instead  of  solving 
for  magnetic  and  electric  fields  directly.  Permittivity  and  permeability  of  free  space 
become  distributed  capacitance  and  inductance,  respectively,  along  the  transmission 
lines. 

Within  the  FD-TLM  method  [1],  space  is  divided  into  a  set  of  three  dimensional 
nodes  as  depicted  in  Fig.  3.1  and  along  the  faces  of  each  cell,  every  E  (electric) 
node  is  surrounded  by  an  H  (magnetic)  field  node  in  any  of  four  directions  and, 
similarily,  every  H  field  node  is  surrounded  by  four  adjacent  E  field  nodes  (Fig  3.2). 
There  are  three  electric  field  nodes  for  each  dimension  of  space  Ex,  Ey,  and  Ez, 
as  well  as  three  magnetic  field  nodes  Hx,  Hy,  and  Hz,  where  each  node  in  the  FD- 
TLM  space  has  a  one-to-one  correspondence  with  points  discretizing  space  containing 
the  structure  to  be  modeled.  The  H-field  and  E-field  nodes  are  connected  through 
transmission  lines  having  characteristic  inductance  and  capacitance  determined  by 
free  space  permittivity  and  permeability.  Each  electric  field  node  connects  the 
two  conductor  transmission  lines  in  a  parallel  fashion  while  the  magnetic  field  nodes 
connect  the  transmission  lines  in  series.  Therefore,  the  E-field  node  is  referred  to  as 
a  shunt  node  (Fig  3.3)  while  the  H-field  nodes  referred  to  as  a  series  node  (Fig  3.4). 

Within  the  FD-TLM  method,  a  variable  mesh  is  used  to  change  the  relative  dis¬ 
tance  between  adjacent  nodes  as  opposed  to  a  uniform  mesh  where  all  of  the  nodes 
are  equidistant  in  space.  This  allows  fewer  cells  to  be  used  for  many  types  of  circuit 
simulations,  and,  therefore  shorter  simulation  time  due  to  fewer  number  of  nodes  to 
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Figure  3.1:  Depiction  of  three-dimensional  space  divided  into  cells. 


Figure  3.2.  Each  cell  in  space  contains  a  set  of  nodes  representing  the  magnetic 
electric  field  nodes.  Note  that  the  nodes  are  on  the  surface  of  the  cube  and 
connected  with  transmission  lines. 


and 

are 
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Figure  3.3:  Depiction  of  the  E-field  shunt  node  where  transmission  lines  are  connected 
in  a  parallel  fashion.  Each  line  represents  two  conductors. 


Figure  3.4:  Depiction  of  the  H-field  series  node  where  transmission  lines  are  connected 
in  a  series  fashion. 
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be  calculated.  To  optimize  simulation  time  and  accuracy,  a  larger  number  of  closely 
spaced  cells  are  used  in  areas  where  fields  are  expected  to  have  a  larger  spatial  gra¬ 
dient,  while  a  fewer  number  of  widely  spaced  cells  are  used  in  the  area,  for  example, 
between  the  circuit  and  the  IC  package.  This  is  analagous  to  an  ordinary  differential 
equations  solver  method  which  uses  smaller  steps  in  the  integration  in  places  where 
the  function  changes  most  rapidly. 

Simulation  within  the  FD-TLM  method  begins  by  defining  a  signal  which  is  ap¬ 
plied  to  an  electric  field  node  or  array  of  electric  field  nodes.  This  signal  propagates 
along  the  transmission  lines  to  the  neighboring  magnetic  field  nodes,  where  part  of 
the  signal  is  reflected  back  to  the  source  E  field  nodes  and  part  is  transmitted  along 
transmission  lines  connected  to  other  E  field  nodes.  The  FD-TLM  method  alternates 
calculating  E  an  H  field  nodes  at  each  time  as  signals  propagating  at  a  cerain  velocity 
travel  from  the  E  field  nodes  to  the  H  field  nodes  and  vice  versa.  The  propagation 
of  signal  pulses  along  the  transmission  lines  and  the  scattering  of  pulses  at  the  nodes 
models  electromagnetic  wave  propagation.  The  voltage  calculated  across  a  shunt  node 
corresponds  to  the  electric  field  at  that  point  in  space  and  the  current  flowing  through 
a  series  node  corresponds  to  the  magnetic  field  at  that  point  in  space. 


3.1  Modeling  Linear  Devices 

Linear  devices  such  as  capacitors,  inductors,  and  resistors  that  maintain  fixed  values 
with  applied  voltage  and  time  are  easily  modeled  in  the  FD-TLM  method  [1],  Note 
the  addition  of  a  permittivity  stub  and  loss  stub  to  the  E-field  node  and  addition  of  a 
permeability  stub  to  the  H-field  node  in  Fig.  3.3,  and  Fig.  3.4.  These  stubs  are  used 
to  describe  conductivity,  permittivity,  and  permeability  at  each  point  in  geometric 
space  and,  therefore,  provide  the  means  for  describing  layers  of  materials,  resistors, 
and  capacitors. 
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Resistors  are  implemented  in  the  FD-TLM  method  by  incorporating  conductivity 
via  the  loss  stubs  for  the  E  field  nodes  in  the  volume  representing  the  resistor.  Con¬ 
ductivity,  <r,  is  calculated  from  resistance,  R,  using  the  relationship  R  =  where 
the  length  of  the  resistor,  L,  is  determined  to  be  in  the  direction  of  current  flow  and 
A  is  the  cross-sectional  area  perpendicular  to  current  flow. 

When  setting  the  conductivity  of  the  nodes  in  space  representing  a  resistor,  the 
region  of  influence  for  each  node  must  be  considered.  That  is,  at  a  point  between 
two  adjacent  nodes,  the  actual  conductivity  is  an  average  of  the  conductivity  at  each 
node.  As  a  result,  while  each  node  in  the  FD-TLM  method  can  be  thought  of  as 
representing  a  cell  with  dimensions  specified  by  the  spacing  between  adjacent  nodes 
in  each  direction  (see  Fig.  3.1).  the  conductivity  within  the  cell  is  not  equal  to  the 
conductivity  of  the  Ex,  Ey,  or  Ez  node,  but,  rather,  is  an  average  of  the  conductivities 
in  adjacent  nodes  or  of  the  nodes  on  the  surface  of  the  cell.  The  actual  value  of 
conductivity  needed  to  implement  a  resistor,  where  (u,v,w)  is  the  dimension  of  the 
(x,  y,  z)  cell,  is  determined  from 


Rx 


4  u 


(TXVW  -f  +  (TX'k—iVWk-.\  +  \ 


(3.1) 


The  subscripts  i,  j,  k  refer  to  the  location  of  each  cell  (see  Fig  3.1)  where  the  H  field 
nodes  on  the  bottom,  left  side,  and  front  of  the  cell  (Fig  3.2)  and  the  E  field  nodes  on 
the  bottom-left,  bottom-front,  and  left-front  of  the  cell  (Fig  3.2)  belong  to  the  (i,j,k) 
cell  and  all  other  nodes  belong  to  adjacent  cells.  The  previous  equation  establishes 
a  resistance  in  the  x  direction;  however,  conductivity  can  be  specified  differently  for 
the  x,  y,  or  z  direction  at  each  cell  allowing  modeling  of  anisotropic  materials. 

Capacitors  are  implemented  by  altering  the  permittivity  at  a  set  of  nodes.  The 
determination  of  permittivity  for  implementation  of  the  capacitor  follows  that  of  the 
resistor.  The  equation  defining  the  capacitance  value  for  an  x-directed  E  field  node  is 
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evw  +  eXtj-iVj- iw  +  eXik-ivvJk-i  +  eXtj-itk-\Vj-iWk-i 

4  u 


(3.2) 


Substrates  and  layers  of  material  are  common  in  ICs,  and,  therefore,  the  FD- 
TLM  method  must  be  capable  of  simulating  these  layers  in  order  to  properly  model 
ICs.  These  layers  are  specified  by  setting  the  permittivity  values  of  the  nodes  in  the 
FD-TLM  method  corresponding  to  the  materials. 


3.2  Modeling  Nonlinear  Devices 

Within  a  non-linear  device,  the  current  and  voltage  do  not  have  a  linear  relationship 
meaning  that  the  conductivity  of  the  device  is  a  function  of  the  applied  voltage. 
Furthermore,  capacitance  for  the  device  may  change  with  applied  voltage  as  occurs  for 
the  diode  or  any  P-N  junction  capacitance.  Due  to  the  relationship  between  voltage 
and  current,  a  non-linear  device  cannot  be  implemented  by  specifying  conductivity 
or  permittivity  in  an  input  data  set  for  the  FD-TLM  method.  Instead,  the  FD-TLM 
method  must  be  given  a  current-voltage  (I-V)  equation  for  the  node  representing 
the  non-linear  device.  The  FD-TLM  method  then  uses  the  I-V  equation  to  update 
the  conductivity  of  the  node  representing  the  non-linear  device  as  a  function  of  the 
voltage  or  electric  field  at  the  node.  Note  that  non-linear  devices  are  implemented  at 
E-field  nodes. 

Many  non-linear  devices,  such  as  the  Josephson  junction,  contain  storage  elements 
which  may  lead  to  a  large  negative  conductivity.  Within  the  FD-TLM  method  a 
voltage  is  calculated  at  each  time  step  and,  from  this  voltage,  the  current  through 
the  non-linear  device  and  its  conductance  is  calculated  and  implemented  at  the  node 
for  the  next  iteration.  However,  an  analysis  of  the  circuit  equivalent  for  a  non-linear 
device  may  reveal  current  sources  which  are  capable  of  driving  current  in  an  active 
convention  rather  than  passive  convention  with  respect  to  voltage  leading  to  negative 
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Figure  3.5:  Equivalent  circuit  of  the  JJ  showing  current  source. 

conductivity  at  the  node.  Within  the  FD-TLM  method,  negative  conductivity  creates 
a  large  instability  and  resultant  oscillations.  As  a  result,  the  storage  elements  of  such 
a  device  cannot  be  simulated  as  a  changing  conductivity  at  the  electric  field  node.  The 
I-V  equation  is  modified,  from  the  circuit  equivalent,  to  contain  all  but  the  storage 
elements  in  question.  The  storage  elements,  specifically  current  sources,  are  modeled 
by  updating  the  electric  field  at  the  electric  field  node  representing  the  nonlinear 
device  to  reflect  the  current.  Figure  3.5,  shows  the  current  source  which  is  modeled 
by  updating  the  value  of  the  E  field  node  at  each  iteration  while  the  resistance 

G(v) 

is  modeled  as  a  conductivity  that  changes  with  time  and  applied  voltage,  ie  the  value 
of  the  loss  stub  is  changed  at  the  E  field  node. 

Non-linear  capacitance  as  well  as  conductance  can  be  implemented  in  the  FD- 
TLM  method  by  changing  the  permittivity  at  the  node  representing  the  device  as  a 
function  of  the  present  value  of  the  electric  field  node.  In  the  case  of  the  Josephson 
junction,  the  capacitance  is  linear  and  is  modeled  as  a  constant  permittivity. 


Chapter  4 

Modeling  the  Josephson  Junction 


4.1  Characteristics  of  the  Josephson  Junction 

Before  beginning  with  implementation  of  the  JJ  and  the  analysis  of  JJ  logic  gates,  it 
is  important  to  understand  the  characteristics  of  the  device.  Figure  4.1  shows  the  dc 
current-voltage  characteristics  of  the  Josephson  junction. 

There  are  two  parts  to  the  dc  curve  which  combine  to  describe  the  superconduct¬ 
ing  tunneling  current  and  single  particle  tunneling  current.  The  superconducting 
tunneling  current  represents  current  in  the  zero  voltage  state  (ZVS)  of  the  device 
where  current  flows  through  the  device  with  no  associated  voltage  drop.  Note  that 
the  current  in  the  ZVS  cannot  exceed  the  critical  current,  Ic,  in  magnitude.  If  the 
current  through  the  junction  should  exceed  the  critical  current  of  the  device,  the  de¬ 
vice  switches  to  another  state  commonly  referred  to  as  the  high  voltage  state  (HVS) 
shown  as  a  jump  from  point  1  to  point  2  in  Fig.  4.1.  A  voltage  develops  across  the 
junction  in  the  HVS  and  follows  the  I-V  characteristics  of  the  single-particle  tunnel¬ 
ing  curve.  Typically,  in  the  HVS,  the  voltage  across  the  junction  will  equal  the  gap 
voltage  V gapi  a  JJ  process  parameters  used  to  describe  the  JJ  current-voltage  rela¬ 
tionship,  with  currents  just  above  the  critical  current  of  the  device.  When  current 
is  decreased  to  a  value  below  the  critical  current,  instead  of  returning  to  the  ZVS, 
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Figure  4.1:  DC  current-voltage  characteristics  of  the  Josephson  junction.  The  vertical 
line  at  the  V=0  line  represents  the  ZVS  and  Cooper  pair  tunneling  while  the  HVS 
curve  or  nonzero  voltage  curve  represents  single-particle  tunneling. 


18 


Figure  4.2:  Equivalent  circuit  for  the  Josephson  junction. 

the  device  remains  in  the  HVS  and  the  I-V  characteristics  continue  to  follow  the 
single-particle  tunneling  curve.  The  device  does  not  return  to  the  ZVS  until  current 
through  the  device  is  lowered  to  zero.  As  a  result,  the  JJ  has  a  hysteresis. 

The  dc  TV  curve  for  the  JJ  is  reconciled  by  the  two  fluid  model  of  supercon¬ 
ductivity  [11].  The  two  fluid  model  provides  the  ability  for  both  superconducting 
Cooper  pairs  and  normal  conducting  electrons  to  contribute  to  current  flow  in  the 
superconductor  at  any  time.  That  is,  normal  conduction  and  superconduction  can 
occur  in  tandem  within  the  superconductor.  This  model  explains  many  phenomena 
within  the  superconductor  including  the  existence  of  resistance  in  the  superconductor 
at  very  high  frequencies,  as  wells  as  how  the  JJ  is  capable  of  switching  between  tun¬ 
neling  of  Cooper  pairs  and  electrons  to  transition  between  a  superconductive  state 
and  a  resistive  state. 

The  circuit  equivalent  for  the  Josephson  junction  is  shown  in  Figure  4.2. 
Insight  into  the  ac  as  well  as  dc  characteristics  of  the  JJ  is  gained  by  examining  the 
equivalent  circuit.  This  consists  of  three  components,  a  current  source  which  repre¬ 
sents  the  tunneling  of  Cooper  pairs,  a  resistance  Qjyj  which  represents  single  particle 
tunneling,  and  a  capacitor  C j  which  is  the  capacitance  created  by  the  tunneling  bar¬ 
rier  sandwiched  between  the  superconductive  electrodes.  As  seen  from  the  equivalent 
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circuit,  superconducting  current  with  its  zero  associated  voltage  drop  can  only  be 
represented  by  the  current  source.  This  is  defined  as 

/  =  Icsin((f>):  (4.1) 

The  JJ  critical  current  is  Ic  and  the  phase  0,  the  quantum  mechanical  phase  difference 
across  the  junction,  is  related  to  voltage  across  the  junction  by 

Tt  =  w  (4'2) 

The  superconducting  current  can  be  rewritten  as 

I  =  Icsin(P0  J  Vdt )  ’  (4.3) 

One  interesting  ac  characteristic  of  the  JJ  can  be  seen  immediately  from  (4.3). 
With  a  constant  applied  voltage,  the  current  through  the  current  source  will  oscillate 
at  a  frequency  a;,  where 

w  <4-4> 

The  JJ  current-voltage  curve  can  be  derived  from  the  equivalent  circuit.  Observing 
Fig.  4.2,  we  start  initially  by  assuming  that  the  current  through  and  the  voltage 
across  the  device  are  zero.  As  an  externally  applied  current  is  increased  through  the 
device,  the  current  first  flows  through  the  capacitor  which  appears  as  a  short.  Then, 
as  a  voltage  is  developed  across  the  capacitor,  current  begins  to  flow  through  the 
resistance  Q^rj-  The  voltage,  generated  by  current  flow  through  the  capacitor  and 
resistance,  changes  the  current  through  the  current  source  by  (4.3).  The  voltage  varies 
as  the  current  begins  to  flow  through  the  current  source  rather  than  the  capacitor 
and  resistor,  and  the  entire  process  continues  until  all  of  the  current  is  absorbed 
by  the  current  source  and  there  is  no  voltage  drop  across  the  junction.  Thus,  the 
JJ  has  remained  in  the  ZVS.  Once  the  system  is  stabilized,  a  change  in  externally 
applied  current  will  repeat  the  entire  process  for  the  new  current  until,  once  again, 
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the  current  source  absorbs  all  of  the  current.  Thus,  a  change  in  externally  applied 
current  generates  a  temporary  spike  in  voltage  until  the  system  stabilizes. 

If  the  externally  applied  current  through  the  junction  is  greater  than  the  critical 
current,  Ic,  of  the  device,  a  voltage  is  generated  across  the  capacitor  and  resistor  as 
before;  however,  the  current  source  cannot  absorb  all  the  externally  applied  current 
lext  in  that  /  =  Icsin(<f> )  can  never  be  equal  to  \ext.  Thus,  some  of  the  externally 
applied  current  must  flow  through  the  resistor,  causing  a  voltage  drop  across  the  JJ. 
The  JJ  is  now  in  the  HVS  and  the  current  source  oscillates  at  the  frequency  given  by 
(4.4).  Once  in  the  HVS,  the  JJ  merely  behaves  as  a  resistor  since  essentially  all  of  the 
externally  applied  current  flows  through  the  resistor.  The  oscillating  current  from  the 
current  source  is  generally  shorted  by  the  capacitor  in  parallel;  however,  some  of  this 
current  shows  up  as  a  small  ripple  in  the  output  voltage  and  current  of  the  device 
when  in  the  HVS. 

Earlier,  it  was  noted  that  a  hysteresis  exists  in  the  dc  characteristics  of  the  JJ.  That 
is,  the  JJ,  once  in  the  HVS,  is  not  capable  of  returning  to  the  ZVS  until  the  current 
through  the  device  is  returned  to  zero.  This  behavior  can  be  explained  considering 
the  ac  characteristics  of  the  JJ,  particularly  that  oscillations  are  created  in  the  current 
source  for  a  constantly  applied  voltage.  When  the  JJ  enters  the  HVS,  a  voltage  of 
approximately  2  mV  is  developed.  This  creates  a  high  frequency  oscillation  in  the 
current  source,  as  mentioned  before,  representing  the  energetic  instability  in  the  JJ 
for  currents  higher  than  the  critical  current  of  the  device.  As  the  externally  applied 
current  decreases  below  the  critical  current,  a  voltage  of  approximately  2  mV  is 
still  observed  across  the  JJ  which  does  not  allow  the  current  source  to  match  the 
externally  applied  current  even  though  the  externally  applied  current  may  be  less 
than  the  critical  current.  The  only  way  to  force  the  JJ  back  into  the  ZVS  is  to 
completely  cut  off  the  oscillations  of  the  current  source  by  reducing  the  externally 
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applied  current  to  zero  and  subsequently  the  voltage  across  the  junction  to  zero. 

4.2  Implementation  of  the  Josephson  Junction  in 
the  FD-TLM  method 

The  equivalent  circuit  of  the  JJ  is  used  to  determine  the  conductance  as  a  function 
of  voltage  to  implementing  the  JJ  in  the  FD-TLM  method.  The  capacitance  for  the 
JJ,  C j  is  not  used  in  the  conductivity  calculation,  but  is  implemented  directly  as  a 
capacitance  within  the  FD-TLM  method  following  the  method  outlined  in  Chapter 
3.  The  following  equations  are  needed  to  determine  the  conductivity  of  the  JJ  as  a 
function  of  voltage  [9] 

I  =  I0sin(<j>)  +  G(V)  +  Cj^-  (4.5) 

PoV  =  ^  (4.6) 

at 

G(V)  =  GiV  +  (/,  +  G,\V\)  |  -  -  -  ,  +  ex‘p  15±r  }  (47) 

where  P0  is  the  plasma  dampening  frequency  of  the  JJ,  <f)  is  the  quantum  mechanical 
phase  difference  for  the  Cooper  pair  particle  across  the  junction  and  changes  with  time 
as  a  function  of  voltage  according  to  (4.2),  I0  is  the  critical  current  of  the  JJ,  and 
G{V)  describes  the  single-particle  tunneling  curve.  A  discussion  of  the  parameters  for 
G(V)  is  obtained  from  [13,  9]  and  are  determined  from  the  measured  properties  of  the 
JJ  after  fabrication.  With  the  current-voltage  relationships  of  the  JJ  determined,  the 
device  can  be  implemented  at  an  E-field  node  or  array  of  E-field  nodes  by  altering  the 
conductivity  at  the  node  or  nodes  as  a  function  of  the  value  of  the  E-field  as  discussed 
in  Chapter  3.  However,  as  was  discovered  during  the  research  stage,  there  are  some 
problems  inherent  in  modeling  the  JJ  as  a  conductance  using  its  equivalent  circuit. 
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An  analysis  of  the  value  of  conductivity  for  the  JJ  with  respect  to  time  reveals 
a  problem  with  simulating  the  JJ  as  just  a  conductivity  which  changes  with  time. 
Initially,  the  FD-TLM  method  begins  by  setting  the  conductivity  of  the  node  repre¬ 
senting  the  JJ  to  zero.  As  signals  develop  within  the  FD-TLM  method,  a  voltage  will 
eventually  be  established  across  the  JJ  node  due  to  the  capacitance  at  the  node  and 
conductivity  can  begin  to  change  accordingly.  The  problem  arises  when  the  JJ  is  in 
the  ZVS  since  current  must  flow  through  the  device  with  no  associated  voltage  drop. 
This,  indeed,  occurs  in  the  FD-TLM  method  as  the  current  source  within  the  circuit 
equivalent  begins  to  absorb  all  of  the  input  current.  However,  the  conductivity  of  the 
node  approaches  infinity  at  this  point.  This  immediately  creates  a  problem  in  any 
numerical  method,  but,  even  if  this  problem  is  corrected  by  imposing  a  maximum 
value  for  conductivity  and  accepting  the  error,  another  problem  still  arises.  With 
an  infinite  conductivity,  any  new  currents  through  the  JJ  will  have  no  effect  on  it 
since  there  is  no  longer  any  way  to  generate  a  voltage  drop  across  the  node,  and,  as 
a  result,  the  JJ  behaves  merely  as  a  shorted  node.  Finally,  it  is  possible  to  have  a 
situation  where  the  current  source  in  the  circuit  equivalent  is  directing  current  in  an 
active  rather  than  passive  convention  with  respect  to  voltage.  In  other  words,  there 
is  energy  storage  in  the  device  and  it  is  possible  for  the  JJ  to  behave  as  a  battery 
leading  to  a  negative  conductivity  at  the  node.  This  creates  instability  within  the 
FD-TLM  method.  Obviously,  a  different  method  must  be  used  for  implementing  the 
JJ. 

Instead  of  simulating  the  entire  device  as  a  node  whose  conductivity  changes  with 
time,  only  the  resistance  is  modeled  as  a  conductivity  at  the  E-field  node.  Ca¬ 
pacitance  is  implemented  as  a  linear  capacitance  using  methods  discussed  in  Chapter 
3. 

The  current  source,  however,  is  incorporated  in  the  FD-TLM  method  following  a 
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procedure  analogous  to  the  implementation  of  a  resonant  tunneling  diode  [12].  That 
is,  modify  the  value  of  the  electric  field  at  the  E-field  node  representing  the  JJ  by 
adding  the  time-averaged  current  source  current  expressed  as 


Is.,,  =  I  +  7s(V”))/ 2, 


(4.8) 


where  Vn  represents  voltage  at  the  present  time  step,  and  Vn+1  is  the  voltage  at 
the  next  time  step  (n  +  l)r  where  r  is  the  FD-TLM  time  step.  Current  Is(Vn)  is 
calculated  as  losing)  where  (f>n  is  the  phase  at  time  nr ,  the  present  time  step,  and  is 
obtained  by  integrating  (4.2)  as  the  sum  of  the  product  P0Vnr  at  each  time  step.  At 
the  start  of  the  simulation,  the  phase,  <f>,  is  initially  zero.  The  new  superconducting 
current,  current  through  the  current  source,  is  found  as 

Is(Vn+1)  =  Is{Vn )  +  ^V+1  ~  Vn).  (4.9) 


Using  a  truncated  Taylor  series,  the  derivative  in  the  above  expression  is 

dIs(V  =  Vn ) 


dV 


=  I0cos((f>n)P0V1' 


yn+ 1  _  ]/n  ■ 


(4.10) 


Incorporating  the  time- average  JJ  superconducting  current  yields  the  new  FD-TLM 
equation  for  updating  the  electric  field  node  voltage 
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+Iosin{cf>n)\ 


(4.11) 


where  V™+1  and  K.n  are  the  shunt  node  voltages  and  the  subscripted  /  are  the  currents 
flowing  in  neighboring  series  (magnetic  field)  nodes  [12].  This  equation  is  for  an 
electric  field  node  for  z-directed  electric  fields.  Analogous  equations  can  be  derived 
for  electric  field  nodes  in  the  x  and  y  directions. 
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Conductivity  for  the  node  is  calculated  as  G{V)  which  will  never  be  infinite  since 
this  represents  normal  current  flow  through  the  resistance  of  the  JJ.  The  current 
source  is  modeled  as  a  current  source  within  the  FD-TLM  method  rather  than  mod¬ 
eling  it  as  part  of  the  conductivity  of  the  electric  field  node  representing  the  JJ. 
Therefore,  once  the  system  is  stable  in  the  ZVS,  all  of  the  current  through  the  JJ  will 
be  represented  by  the  current  source  while  maintaining  a  reasonable  conductivity  of 
the  node  as  being  (?(0)  which  is  just  the  Rpap  resistance,  where  R5ap  is  the  equivalent 
resistance  of  the  JJ  in  the  HVS.  For  any  additional  applied  current,  the  FD-TLM 
method  responds  by  creating  voltages  across  the  capacitance  and  G(V )  to  readjust 
the  current  source  for  the  newly  applied  current.  To  summarize,  this  JJ  model  imple¬ 
mentation  never  creates  a  short  at  the  E  field  node  and  thus  never  becomes  trapped 
in  the  ZVS,  and  lacks  the  instability  of  having  negative  conductivity,  solving  all  of 
the  problems  associated  with  the  original  flawed  modeling  approach. 

Figure  4.3  shows  the  typical  simulation  of  a  single  JJ  device  with  the  FD-TLM 
method.  In  the  simulation,  a  5  mV  pulse  is  in  series  with  a  27.0  ohm  resistor  and  the 
JJ,  providing  current  greater  than  the  critical  current  of  the  JJ.  The  voltage  level  and 
current  through  the  circuit  rises  until  the  JJ  enters  the  HVS  and  an  output  voltage 
is  appears  across  the  junction.  Furthermore,  the  externally  applied  current  decreases 
once  the  JJ  is  in  the  HVS  showing  that  the  JJ  has  switched  from  a  superconducting 
state  to  a  resistive  state.  The  externally  applied  current  is  created  by  a  voltage 
source  in  series  with  a  resistor.  Thus,  when  there  is  a  non-zero  voltage  across  the  JJ, 
the  circuit  current  decreases  according  to  Kirchhoff’s  voltage  law.  The  characteristic 
ripple  created  by  the  oscillating  current  source  in  the  HVS  can  also  be  observed  in 
both  the  voltage  across  the  JJ  and  in  the  current  through  the  JJ.  Next,  the  voltage 
source  is  cut  off  in  the  circuit  reducing  the  current  to  zero  in  the  circuit  while  the  JJ 
requires  time  to  leave  the  HVS  and  settle  back  into  the  ZVS.  Comparison  to  references 


Figure  4.4:  Depiction  of  the  arrangement  of  JJs  within  the  SQUID  loop. 

[13,  9],  show  that  this  is,  indeed,  the  correct  general  behavior  for  the  JJ.  The  single 
JJ  was  also  modeled  using  conventional  circuit  simulation  performed  by  solving  the 
differential  equations  for  the  JJ  using  a  fifth-order  Runge-Kutta  method  (use  of  this 
method  is  discussed  in  the  next  chapter).  Fig.  4.3  reveals  good  agreement  between 
FD-TLM  and  conventional  simulation  validating  the  FD-TLM  JJ  modeling  approach, 
results  of  conventional  circuit  simulation 

With  the  implementation  of  the  JJ  devices  discussed  so  far,  a  wide  variety  of  JJ 
logic  circuits  can  be  simulated  with  the  FD-TLM  method.  However,  when  JJs  are 
placed  within  a  superconducting  loop  as  depicted  in  Fig.  4.4,  the  two  JJs  become 
coupled  as  a  result  of  the  requirement  that  the  flux  through  a  superconducting  loop 
is  quantized.  This  configuration  is  referred  to  as  a  Superconducting  Quantum 
Interference  Device  (SQUID).  Within  the  SQUID  loop,  the  JJs  are  coupled  via  the 
vector  magnetic  potential  existing  within  the  superconducting  material  forming  the 
loop  connecting  the  two  JJs,  or  equivalently,  via  the  magnetic  flux  passing  through 
the  loop.  This  flux  forces  a  fixed  phase  change  difference  between  the  JJs  to  quantize 
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the  flux  through  the  loop.  Therefore,  instead  of  behaving  independently,  the  phase  of 
one  JJ  becomes  a  function  of  the  phase  across  the  other  JJ.  Since  the  phase  across  the 
junction  determines  the  current  through  the  junction  as  given  by  (4.1),  the  current 
flowing  through  one  JJ  becomes  a  function  of  the  current  traveling  through  the  other 
JJ  and  is  also  a  function  of  the  magnetic  flux  passing  through  the  SQUID  loop. 

Quantum  mechanics  must  be  used  to  determine  the  relationship  between  the  JJ 
currents  and  the  applied  external  magnetic  flux.  Consider  a  Cooper  pair  wave  function 
that  circulates  through  the  SQUID  loop  and  then  integrate  the  phase  of  this  wave 
function  around  the  loop  as  described  in  [13].  This  yields  the  following  equation 
where  a  discontinuity  of  the  phase  occurs  at  the  junction  interface  yielding, 

j>  \79  ■  dl  =  2mr  =  9\  +  62  (4-12) 

where  6\  is  the  phase  across  Ji  and  62  is  the  phase  across  J2,  and  the  integration  of 
the  phase  around  the  loop  must  be  an  integer  multiple  of  2n.  Through  the  use  of  the 
gauge-invariant  phase  we  obtain 

2e  r 

O2  =  +  —  j>  A  •  dl  —  2mr,  (4-13) 


which  gives  a  relationship  between  phase  and  the  magnetic  potential  vector.  Use  of 
Stokes’  theorem  relating  the  magnetic  vector  potential  to  flux  gives  the  final  equation 
for  expressing  the  phase  relationship  between  Ji  and  J2  for  an  externally  applied  flux 

9  7 r  (f) 

*2  =  *1  -  (4-14) 

where  $e  is  the  externally  applied  flux  through  the  loop  and  $0  is  the  flux  quan¬ 
tum  and  is  equal  to  2.068xl0-15  Wb/m2.  The  total  current  through  the  SQUID  is 
obtained  from  the  following  expression 


It  =  Icisin(9i)  +  IC2sin(91  - 


2?r$e 

$0  h 


(4.15) 
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Figure  4.5:  Total  critical  current  versus  applied  flux  for  the  symmetric  SQUID. 

which  is  used  to  determine  the  total  critical  current  of  the  SQUID  as  a  function  of 
the  applied  flux  by  finding  the  phase  which  maximizes  the  total  current  at  each  value 
of  flux.  For  the  case  of  the  SQUID  where  both  JJs  are  identical  and  the  SQUID  loop 
is  symmetric,  a  simple  expression  can  be  obtained  for  the  total  critical  current 


lTc($e)  =  2/ci 


7 r$e 
cos - 


(4.16) 


A  plot  of  the  total  critical  current  versus  applied  flux  is  given  in  Fig.  4.5,  in  which  the 
interference  pattern  formed  classifies  the  SQUID  as  an  interference  device.  Note 
that  the  critical  current  versus  applied  flux  oscillates  with  increased  flux  rather  than 
changing  linearly.  An  applied  flux  equal  to  one-half  of  the  flux  quantum  yields  a  total 
critical  current  of  zero,  and  a  flux  equal  to  an  integer  multiple  of  the  flux  quantum 
yields  the  maximum  critical  current.  Physically,  this  occurs  due  to  the  adjustment  of 
•J  J  currents  to  keep  an  integer  multiple  of  the  flux  quantum  flowing  through  the  loop. 

An  analytical  expression  for  total  critical  current  versus  externally  applied  flux 
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was  obtained  for  the  case  of  the  symmetric  loop.  A  symmetric  loop  is  one  where  the 
Josephson  junctions  have  the  same  properties  and  are  placed  symmetrically  within 
the  loop.  However,  in  the  case  where  the  loop  is  assymmetric,  it  is  difficult  to  derive 
an  analytical  expression  for  total  critical  current.  Instead,  a  numerical  technique  must 
be  used  which  directly  solves  for  the  maximum  allowable  current  through  the  SQUID 
at  each  value  of  applied  magnetic  flux. 

It  has  been  assumed  that  currents  through  the  loop  do  not  contribute  to  the 
magnetic  flux  through  the  loop  when,  in  fact,  they  do.  Inclusion  of  the  contribution 
to  magnetic  flux  in  the  loop  from  currents  in  the  loop  leads  to  a  new  expression  for 
the  JJ  phases  versus  applied  flux 

e2  =  «,  -  +  (4.i7) 

^  O 

where  Li  is  the  self-inductance  for  the  side  of  the  SQUID  loop  containing  Jx,  L2  is 
the  self-inductance  for  the  side  of  the  loop  containing  J2,  U  is  the  current  through 
Li,  and  I2  is  the  current  through  L2  [13]. 

To  make  the  expression  for  total  critical  current  of  the  SQUID  more  complete, 
the  effect  of  finite  junction  size  can  be  included.  Previously,  we  have  assumed  for 
purposes  other  than  the  calculation  of  critical  current  density  and  capacitance,  that 
the  junction  is  a  point  junction.  However,  with  finite  area,  it  is  possible  for  magnetic 
flux  to  exist  in  the  JJ  tunneling  barrier,  altering  the  critical  current  of  the  JJ  via  the 
expression 

(4.18) 

where  $  is  the  flux  contained  within  the  Josephson  tunneling  area. 

Considering  these  nonidealities,  the  expression  for  determining  the  total  critical 
current  for  the  SQUID  loop  as  a  function  of  externally  applied  magnetic  flux  can 
become  unwieldy,  requiring  a  numerical  method  for  solution.  This  is  where  the  FD- 
TLM  method  becomes  advantageous. 


Ic{$)  =  /c(0) 


sin  (n  */*o) 


(7T$/$0) 
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To  implement  the  SQUID  in  the  FD-TLM  method,  a  new  model  is  created  which 
will  couple  two  JJs  using  the  methods  just  described.  That  is,  one  JJ  within  the 
SQUID  loop  is  allowed  to  vary  freely  as  if  acting  alone  according  to  the  relationships 
(4.5, 4. 6, 4. 7).  However,  the  second  JJ  in  the  loop  specified  with  SQUID  parameters 
given  to  the  FD-TLM  method  adjusts  its  phase  based  upon  a  calculation  of  the 
flux  through  the  surface  area  inside  the  loop  and  the  phase  of  the  JJ  to  which  it 
is  coupled.  Here,  instead  of  determining  a  function  of  total  critical  current  versus 
applied  flux  and  then  treating  the  SQUID  as  a  single  JJ  as  is  done  in  many  methods 
[9],  a  coupling  relation  between  the  JJs  within  the  loop  is  directly  implemented.  Since 
the  value  for  externally  applied  flux  is  calculated  by  the  summation  of  H  field  nodes, 
the  self-induced  flux  described  earlier  is  included  since  currents  through  the  loop 
in  the  FD-TLM  method  will  automatically  create  magnetic  fields  as  a  result  of  the 
solution  of  Maxwell’s  equations. 


Chapter  5 

Resistively  Coupled  Josephson 
Logic  Circuits 

5.1  Introduction 

In  the  previous  chapter,  all  of  the  concepts  needed  to  implement  a  single  JJ  or  a 
SQUID  were  covered  along  with  a  method  for  implementing  the  devices  within  the 
FD-TLM  method.  In  this  and  the  next  chapter,  emphasis  is  placed  on  logic  circuits 
utilizing  the  JJs  and  validation  of  the  simulation  results  from  the  FD-TLM  method. 
The  FD-TLM  method  as  discussed  in  Chapter  3  solves  Maxwell’s  equations  for  elec¬ 
tric  and  magnetic  fields  points  in  discretized  three-dimensional  space  at  each  time 
step  within  the  numerical  method.  Within  the  discretized  three-dimensional  space, 
properties  such  as  conductivity,  permittivity,  and  permeability  are  specified  for  each 
point.  As  a  result,  to  simulate  a  circuit,  the  FD-TLM  method  requires  the  physical 
properties  and  dimensions  of  the  circuit  as  it  would  be  fabricated.  A  CAD  system 
such  as  MAGIC  [4]  could  be  altered  to  create  a  data  set  necessary  for  the  FD-TLM 
method  allowing  a  chip  design  to  be  analyzed  quickly  without  component  parasitic 
extraction. 

One  technique  for  simulating  JJ  logic  circuits  is  to  perform  a  SPICE  [3]  analysis 
with  component  values  extracted  from  the  physical  layout  of  the  circuit.  Since  a 
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Figure  5.1:  The  JAWS  logic  gate. 

version  of  SPICE  modified  to  model  JJs  was  not  readily  available,  a  conventional 
circuit  analysis  based  on  Kirchkoffs’  current  law  is  used  instead.  It  will  be  found  that 
the  FD-TLM  method  not  only  obtains  good  results,  but  often  predicts  the  existence 
of  additional  parasitics  or  other  phenomena  omitted  in  conventional  circuit  analysis. 


5.2  The  JAWS  circuit 

The  Josephson-Atto  Weber  Switch  (JAWS)  [14]will  be  the  first  JJ  logic  circuit  to 
be  analyzed  with  the  FD-TLM  method.  This  circuit  is  commonly  referred  to  as 
resistively  coupled  Josephson  logic  (RCJL)  in  that  the  circuit  switches  output  states 
from  “low”  to  “high”  by  increasing  the  amount  of  current  injected  into  the  JJ  devices 
beyond  their  critical  current.  A  better  understanding  is  obtained  by  analyzing  the 
circuit  for  the  JAWS  logic  gate  (Fig.  5.1).  Initially,  both  JJs  within  the  circuit 
are  in  the  ZVS  and  appear  as  “shorts”.  A  steady-state  biasing  current  is  established 
by  Vi0  and  Rso  in  the  circuit  which  flows  through  J2  to  ground.  Voltage  inputs  Va 
and  Vf,  are  applied  at  the  input  resistances  Ra  and  R&,  respectively.  Voltage  applied 
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to  the  inputs  creates  additional  current  through  J2.  If  the  sum  of  the  input  currents 
and  bias  current  exceeds  the  critical  current  of  J2,  J2  will  proceed  to  the  HVS  where 
it  behaves  as  a  resistance  rather  than  a  short.  At  this  point,  the  lfi  resistor,  R^,  will 
have  a  much  smaller  resistance  than  the  equivalent  resistance  of  J2  and  the  biasing 
current  will  flow  through  Ji  and  then  to  ground  through  Rj  while  the  input  currents 
flow  through  R^.  Note  that  although  J2  is  in  the  HVS  at  this  point,  there  is  still  no 
output  current  flowing  through  R0Uf.  This  leads  to  the  purpose  of  Ji.  Ji  has  a  critical 
current  designed  so  that  the  sum  of  the  input  currents  flowing  through  Ji  will  not 
force  Ji  into  the  HVS,  while,  in  contrast,  the  biasing  current  flowing  through  Ji  will 
force  Ji  into  the  HVS.  Therefore,  with  the  biasing  current  flowing  through  Ji  after  J2 
has  reached  the  HVS,  Ji  proceeds  to  the  HVS,  forcing  the  input  current  to  continue 
flowing  through  R^  while  the  biasing  current  now  flows  through  R^.  As  a  result, 
the  JJ  logic  circuit  or  JAWS  circuit  has  gated  “high”.  To  return  to  a  logic  “low”,  the 
bias  current  and  input  currents  must  all  be  returned  to  zero  due  to  the  hysteresis  of 
the  JJ. 

Josephson  junction  logic  circuits  operate  on  the  principle  of  switching  current 
among  different  paths.  The  JAWS  gate  is  considered  a  resistively  coupled  logic  gate 
because  it  relies  on  forcing  a  current  through  the  JJ  greater  than  its  critical  current 
to  switch  logic  levels.  Later,  a  discussion  of  magnetically  coupled  logic  (MCJL)  gates 
will  clarify  the  distinction  between  RCJL  and  MCJL.  Nevertheless,  in  this  circuit, 
J2  performs  all  the  work  in  terms  of  switching  the  logic  gate  while  Ji  isolates  the 
input  from  the  output,  forcing  the  input  current  through  R<*  and  the  biasing  current 
through  the  output  resistor. 

The  JAWS  logic  gate  can  operate  as  either  an  AND  gate  or  an  OR  gate  depending 
on  the  value  of  the  biasing  current  established  through  J2.  This  is  performed  by 
establishing  the  bias  current  using  the  relationship  Ia  +  h  +  hias  >  hi  where  Ic  is 
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the  critical  current  of  J2,  and  Ia  and  If,  are  the  input  currents.  The  gate  behaves  as 
an  OR  gate  when  Ia  +  hias  >  h  or  h  +  hias  >  h  and  behaves  as  an  AND  gate  only 
when  Ia  +  h  +  hias  >  h  but  with  the  constraints  h  +  hias  <  h  and  Ia  +  hias  <  h- 

5.2.1  FD-TLM  implementation 

Th e  fdtgraph  [16]  program  was  developed  to  create  data  sets  for  simulation  of  the  var¬ 
ious  JJ  logic  circuits  in  the  FD-TLM  method.  Fdtgraph  is  a  mouse-driven  graphical 
interface  allowing  the  user  to  graphically  layout  the  circuit.  Fdtgraph  then  automat¬ 
ically  creates  the  ASCII  text  data  set  required  for  the  FD-TLM  simulation  of  the 
circuit  layout. 

The  JAWS  circuit  is  simulated  within  a  50x50x50  /am3  (x,y,z)  box  filled  with  air 
dielectric  with  a  grid  spacing  (spacing  between  points  in  space)  of  1  /im.  The  walls 
of  the  box  are  modeled  using  perfect  conductors  as  boundary  conditions,  and  the 
bottom  of  the  box  is  used  as  a  superconducting  ground  plane. 

Fig.  5.2  shows  the  typical  fabrication  steps  for  JJ  integrated  circuits.  This 
structure  is  simulated  within  the  FD-TLM  method;  however,  since  the  silicon  sub¬ 
strate  is  used  as  a  base  for  construction  and  electrical  properties  are  the  primary 
concern  within  circuit  simulation,  the  silicon  substrate  is  excluded  from  the  layout 
within  the  FD-TLM  method.  Thus,  the  bottom  surface  of  the  box  containing  the 
circuit  is  used  as  the  superconducting  ground  plane  and  a  2  ym  thick  layer  of  Si02 
is  used  as  the  insulator  between  the  ground  plane  and  base  electrodes  or  first  layer 
of  metal  for  the  JJs.  The  tunneling  barrier  is  implemented  as  a  y-directed,  the 
field  node,  which  has  current-voltage  characteristics  representing  that  of  the  single 
JJ,  centered  in  the  area  representing  the  JJ  tunneling  material  as  shown  in  Fig.  5.3. 

Superconducting  wiring  is  implemented  as  infinitely  thin,  perfectly-conducting  seg¬ 
ments  of  metal  within  the  FD-TLM  method.  Interconnects  are  made  intentionally 
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Figure  5.2:  Cross-section  of  layers  used  in  fabrication  of  JJ  ICs. 

I 

’  i 


Ey  node  implementing  the  Josephson  junction 


Figure  5.3:  Depiction  of  the  actual  placement  of  the  Ey  node  in  the  tunneling  barrier 
used  to  simulate  the  JJ. 
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Figure  5.4:  JAWS  layout  used  for  FD-TLM  simulation. 

short  to  reduce  reflection  and  electromagnetic  noise  enabling  comparison  with  conven¬ 
tional  circuit  analysis.  The  FD-TLM  layout  is  shown  in  Fig.  5.4.  The  parameters 
are  implemented  following  [12]  Ic  =  0.1  mA,  Gt  =  4.0  mS,  h  =0  A,  G2  =0.1  S, 
V.  =2.0  mV,  Vt  =0.1  mV,  P0  =3.039  x  1015  Wb"1,  and  Cj  =0.5  pF  for  the  JJ.  Three 
voltage  sources  are  implemented  as  shown  in  Fig.  5.1. 

The  entire  simulation  is  performed  for  350  ps  with  a  1.7  ps  time  step.  The  results 
of  the  FD-TLM  simulation  are  shown  in  Fig.  5.5. 

The  biasing  current  was  established  to  allow  the  JAWS  logic  gate  to  behave  as 
an  AND  gate.  In  Fig.  5.5,  the  biasing  current  is  established  by  V50  and,  as  soon  as 
the  system  is  stable,  an  input  voltage  of  2  mV  is  applied  to  input  A.  At  this  point, 
the  output  voltage  remains  at  zero  as  J2  absorbs  all  of  the  input  current  without 
triggering  to  the  HVS.  At  a  later  time,  while  leaving  the  input  at  A  “high”  a  2  mV 
voltage  is  also  applied  to  input  B.  At  this  point,  with  some  delay,  J2  is  triggered  into 
the  HVS  and  the  gate  goes  “high”.  A  signal  of  approximately  1.7  mV  is  observed  at 
the  output  which  can  be  used  as  the  input  for  a  subsequent  gate.  From  the  results 


of  simulation,  it  is  seen  that  the  FD-TLM  simulation  method  provided  the  correct 
general  behavior  of  the  circuit. 


5.2.2  Validation  of  FD-TLM  results 

Conventional  circuit  analysis  is  used  to  validate  the  FD-TLM  results  and,  for  the 
Josephson  junction  circuit,  involves  a  three  step  process  noting  that  SPICE  version 
3f.4  cannot  model  a  JJ.  Nodal  analysis  is  performed  on  the  circuit  followed  by  a  run 
through  MAPLE  to  solve  for  a  set  of  first  order  differential  equations  which  are,  then, 
solved  using  an  numerical  ordinary  differential  equation  integration  method. 

First,  nodal  analysis  is  performed  on  the  circuit  to  obtain  a  set  of  differential 
equations  governing  operation  of  the  circuit.  The  following  are  the  nodal  equations 
for  the  JAWS  circuit. 


^  +  ^  +  ^  +  £  +  /i  =  ° 

-h  +  C3d-^  +  G\ ( V2  -  V3)  +  Iasinifa)  =  0 
^  +  CjM  +  Gl(V3)  +  IC2sin(4>2)  -  C3d-^ 

—Gi(V2  —  Vs)  —  Icisin(<f>i)  +  I2  =  0 
~h  +  CoUT^t  +  =  0 

djt  =  Po(V2  -  V3) 

dj£  =  Po(v3) 

Vs-V4  =  L2% 

MAPLE  [4]  is  then  used  to  solve  for  the  first  order  derivatives  of  each  variable  (see 
Fig.  5.6)  from  (5.1). 

Using  Borland  C/C++  version  4.0  on  a  PC,  a  program  was  written  using  the 
fifth-order  Runge-Kutta  adaptive  time  step  routine  [5]  to  solve  the  set  of  first-order 
differential  equations  (Appendix  B).  The  fifth-order  Runge-Kutta  method  is  the  pre- 
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>  ans[l]; 

|phil(0  =  /MV2(/)-V3(0) 

>  ans  [  2  ] ;  ~  ~  ”  ~ .  “  . 

|phi2(/)  =  ^V3(0 

>  ans[3] ;  ~ ~ 

iv,M.  2 RdVl(t) - RdVt(t) - RdVbjt)  +  VI (t) Rin  +  11(f) Rin Rd 
st  H;  : 

>  ans[ A ] ;  ~ 

0 

“  V2(/)  =  -(-2  Il(/)/&o  +  V3(/)-  Vso(/)  +  Gl(V3(/))7iso  +  /c2sin(phi2(0)^ 

+ 12 (/)  /fao  +  Gl(  V2(0  -  V3(r))  +  Id  sin( phi !(/))  Rso)/(Cj  Rso) 

>  ans  [  5  ] ;  —  “  ~  ™“~  ‘ ' 

d  v2((\-  jUO_^l±  V3(Q-  Vso(0  +  Gl(V3(/))/fro  +  /c2sin(phi2(Q)/foo  +  I2(r)flso 
V  ;  Cj Rso 

>  ans  [  6  ] ;  *  ~  ~  "  ~  "  ’  ~~~  ~ ~  “ 


>  ans[7] 

>  ans [8] 


8  VVn_V(t)RQut-V4(t) 
St  Cout  Rout 


-Vl«)  +  V2«) 
LI 


V3(Q-V4(Q 

L2 


Figure  5.6:  MAPLE  is  used  to  solve  first  order  derivatives. 
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ferred  method  for  solving  the  ordinary  differential  equations,  in  this  instance,  since 
lookup  tables  are  used  to  determine  the  values  of  input  voltages  at  each  time  step  and 
the  input  voltage  curves  are  not  smooth.  The  results  for  the  Runge-Kutta  method 
applied  to  the  characteristic  differential  equations  of  the  JAWs  circuit  are  shown  in 
Fig.  5.5  as  well  as  the  comparison  to  the  simulation  of  the  JAWS  circuit  with  the 
FD-TLM  method. 

During  the  process  of  simulating  the  JAWS  circuit  with  the  Runge-Kutta  method, 
it  was  discovered  that  a  few  parasitics  automatically  modeled  by  the  FD-TLM  method 
were  not  included  in  the  conventional  circuit  analysis.  Upon  analysis  of  the  differing 
results  between  the  two  methods,  the  inclusion  of  two  inductors  Li  and  L2  within  the 
conventional  circuit  analysis  resulted  in  very  close  agreeement  between  the  FD-TLM 
and  Runge-Kutta  methods.  The  physical  dimensions  for  the  segments  of  conductors 
important  for  parasitic  inductance  determination  were  fed  to  the  FastHenry  program 
[2]  which  gave  realistic  values  for  inductances  listed  in  Fig.  5.1.  Capacitances  C d  and 
Cout  were  added  to  the  circuit  to  represent  capacitance  between  the  circuit  and  ground 
plane  through  the  Si02  insulating  layer.  In  this  case,  it  was  discovered  that  the  FD- 
TLM  method  often  predicts  the  existence  of  parasitics  which  might  be  overlooked  in 
conventional  circuit  analysis. 

Slight  differences  between  the  FD-TLM  and  conventional  circuit  analyses  method 
can  be  attributed  to  slightly  inaccurate  extraction  of  the  parasitic  component  values. 
For  example,  within  the  FastHenry  inductance  calculation  program,  the  conductors 
were  assumed  to  be  0.1  microns  thick  while  the  conductors  are  modeled  as  being 
infinitely  thin  in  the  FD-TLM  method.  Furthermore,  instead  of  describing  the  entire 
superconductor  wiring  layout,  only  the  physical  dimensions  and  locations  for  conduc¬ 
tors  considered  most  important  in  determining  parasitic  inductance  were  included 
in  the  FastHenry  data  set.  There  is  very  close  agreement  between  the  FD-TLM 


and  circuit  analyses,  and  the  results  agree  with  the  general  behavior  expected  which 
validates  the  FD-TLM  method  approach  for  simulating  RCJL  logic  circuits. 


Chapter  6 

Magnetically  Coupled  Josephson 
Logic  Circuits 

6.1  Introduction 

In  chapter  5,  FD-TLM  simulations  of  resistively  coupled  Josephson  logic  circuits  were 
performed  and  verified.  However,  as  discussed  in  Chapter  4,  MCJL  circuits  operate 
on  a  different  principle  than  RCJL  circuits  requiring  a  different  simulation  method  to 
include  the  dependence  of  critical  current  for  the  J  Js  as  a  function  of  externally  applied 
flux.  This  chapter  will  show  that  the  FD-TLM  method  is  capable  of  implementing 
magnetically  coupled  Josephson  logic  circuits.  Two  different  types  of  logic  circuits 
will  be  analyzed,  the  DC  SQUID  and  the  (Modified  Variable  Threshold  Logic)  MVTL 
[14]  circuit. 

6.2  DC  SQUID  Simulation 

The  circuit  for  the  2  JJ  DC  SQUID  is  shown  in  Fig  6.1.  This  circuit  is  referred  to 
as  magnetically  coupled  Josephson  logic  (MCJL)  as  opposed  to  resistively  coupled 
Josephson  logic  (RCJL)  discussed  for  the  JAWS  circuit.  The  difference  is  that  the 
MCJL  circuit  uses  magnetic  coupling  to  switch  between  the  HVS  and  ZVS  instead 
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Figure  6.1:  Equivalent  circuit  for  the  2  JJ  DC  SQUID. 

of  current  injection.  This  will  be  discussed  in  detail.  The  fundamental  part  of 
the  two  JJ  DC  SQUID  logic  gate  is  the  SQUID  loop  which  operates  as  discussed  in 
Chapter  3.  In  Fig.  6.1,  a  biasing  current  is  established  by  Vso  with  resistor  Rso  which 
will  flow  through  the  SQUID  loop  since,  initially,  both  JJs  are  in  the  ZVS  and  appear 
as  superconducting  shorts. 

Later,  an  input  voltage  Vctn  is  applied  to  the  control  line  creating  a  current  and, 
thus,  a  magnetic  flux.  The  control  line  is  magnetically  coupled  to  the  SQUID  and 
creates  the  external”  flux  which  changes  the  total  critical  current  of  the  SQUID  loop. 
Switching  of  the  SQUID  to  the  HVS  follows  the  basic  procedure  discussed  in  Chapter 
3.  If  the  current  through  the  control  line  and,  therefore,  the  externally  applied  flux, 
is  high  enough,  the  SQUID  will  no  longer  be  capable  of  supporting  the  bias  current 
as  a  superconducting  current  and  the  JJs  within  the  loop  will  change  to  the  HVS. 
As  a  result,  the  bias  current  will  be  shunted  through  the  output  resistance  R 0ut  and 
the  gate  will  be  in  the  “high”  voltage  state.  To  reset  the  SQUID  to  the  ZVS  (logic 
“low”)  from  the  HVS,  the  bias  current  established  by  Vso  must  be  reduced  to  zero. 
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Figure  6.2:  Physical  layout  of  the  DC  SQUID. 

6.2.1  FD-TLM  implementation 

Once  again,  fdtgraph  is  used  to  create  the  layout  of  the  circuit  and  the  data  set  for  the 
FD-TLM  method.  The  circuit  layout  is  seen  in  Fig.  6.2.  A  50x50x50  ^m3  (x,y,z) 
metal  box  is  filled  with  air  dielectric  with  perfectly-conducting  boundary  conditions 
for  the  walls  of  the  box.  Using  the  typical  fabrication  technique  for  JJ  circuits 
the  bottom  of  the  box  is  used  as  the  superconducting  ground  plane  while  a 
thick  layer  of  SiCU  is  used  between  the  ground  plane  and  the  base  electrodes.  All 
of  the  conductors  are  modeled  as  1  /urn  thick  perfect  conductors  and  are  generally  5 
fim  wide.  The  JJ  parameters  are  the  same  as  for  the  JAWS  circuit  discussed  in  the 
previous  chapter  and  each  JJ  is  implemented  as  a  single  node  centered  in  the  area 
to  representing  the  tunneling  barrier  as  before.  Resistors  are  implemented  following 
that  of  Chapter  3.  Finally,  voltage  sources  are  implemented  as  shown  in  Fig.  6.1. 
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Figure  6.3:  Results  of  simulation  of  the  DC  SQUID  logic  gate. 

The  simulation  duration  was  185ps  and  the  time  step  was  1.7  fs.  The  results  of  the 
simulation  are  shown  in  Fig.  6.3. 

During  the  simulation,  as  shown  in  Fig.  6.3,  a  voltage  signal  Vso  is  applied  to 
generate  a  bias  current  through  the  SQUID  loop.  At  this  time,  the  output  voltage 
remains  zero,  and,  at  a  later  time,  a  voltage  is  applied  to  the  control  line,  ~Vcin  creating 
a  current.  With  a  small  delay,  the  output  voltage  increases  to  approximately  1.7  mV. 
When  the  input  signal,  Vcin  returns  to  zero,  the  output  voltage  remains  at  1.7  mV 
until  Vso  returns  to  zero.  This  follows  the  expected  general  behavior  of  the  circuit 
validating  the  FD-TLM  method. 

6.2.2  Validation  of  FD-TLM  results 

The  same  procedure  used  for  the  JAWS  circuit  is  used  to  validate  the  FD-TLM  results 
for  the  two  JJ  DC  SQUID.  The  circuit  with  labeled  nodes  as  shown  in  Fig.  6.1  is 
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used  for  nodal  analysis  and  to  obtain  the  following  set  of  equations. 
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All  inductance  values  and  mutual  inductances  were  extracted  from  the  layout  by 
providing  a  description  of  the  layout  to  the  FastHenry  program  (see  Appendix  C). 
MAPLE  was  then  used  to  solve  the  nodal  equations  for  a  set  of  first-order  differential 
equations  (see  Appendix  C).  The  fifth-order  adaptive-time  step  Runge-Kutta  method 
was  used  to  solve  the  set  of  differential  equations,  the  results  of  which  are  shown  in 
Fig.  6.3.  along  with  comparison  to  the  results  of  the  FD-TLM  method.  Excellent 
agreement  between  the  two  methods  validates  the  FD-TLM  method. 

Several  conventional  circuit  simulations  of  the  SQUID  were  performed  to  judge 
the  sensitivity  of  the  circuit  performance  to  changes  in  the  values  of  parasitic  compo¬ 
nents.  Capacitance  to  ground  and  capacitance  between  the  control  line  and  SQUID 
loop  were  included  in  prior  simulations  and  were  found  to  have  little  effect  on  cir- 
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Figure  6.4:  Equivalent  Circuit  for  the  MVTL  logic  gate. 

cuit  operation  while  different  inductances  values  greatly  altered  the  operation  of  the 
circuit.  Minor  differences  in  the  results  for  the  FD-TLM  method  and  the  conven¬ 
tional  circuit  simulation  method  can  be  attributed  to  slightly  inaccurate  values  for 
the  inductances  extracted  from  the  physical  layout  of  the  circuit. 

6.3  MVTL  Simulation 

The  MVTL  circuit  uses  a  combination  of  current  injection,  increasing  current 
through  the  JJ  to  trigger  the  onset  of  the  HVS,  and  magnetic  coupling  to  lower  the 
effective  critical  current  of  the  SQUID  loop.  Fig.  6.4  shows  this  clearly,  where  the 
input  current  flows  through  the  control  line  creating  flux  to  lower  the  total  critical 
current  of  the  SQUID  and,  at  the  same  time,  is  injected  into  the  the  top  of  the 
SQUID  loop  to  increase  the  current  through  the  SQUID  loop  and  force  a  more  rapid 
transition  to  the  HVS.  The  combination  of  current  injection  and  magnetic  coupling 
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is  the  key  to  giving  the  MVTL  logic  gate  the  fastest  switching  time  of  all  the  logic 
gates  discussed. 

Analyzing  the  operation  of  the  MVTL  circuit,  a  biasing  current  through  the 
SQUID  loop  is  established  by  the  voltage  source  Vso.  At  this  point,  the  total  current 
through  the  SQUID  loop  is  less  than  the  total  critical  current,  and  the  JJs  remain 
in  the  ZVS,  behaving  as  superconducting  shorts.  When  an  input  signal  is  applied  to 
the  control  line,  the  current  creates  a  flux  which,  when  coupled  to  the  SQUID  loop, 
lowers  the  total  critical  current  of  the  SQUID  loop.  The  input  current  flows  to  ground 
through  J3  and  the  SQUID  loop,  which  are  both  in  the  ZVS.  The  reduction  of  total 
critical  current  and  injection  of  additional  current  sends  the  SQUID  into  the  HVS. 
At  this  point,  the  biasing  current  will  flow  through  J3  and  Rccmj  to  ground,  and  the 
input  current  caused  by  Vcin  will  flow  through  Rcou(  to  ground.  The  biasing  current 
has  a  larger  value  than  the  input  current  and  will  force  J3  into  the  HVS,  which  then 
steers  the  input  current  through  Rcout,  while  the  biasing  current  flows  through  the 
output  resistance  R out.  At  this  point,  the  MVTL  circuit  is  in  the  logic  “high”  state. 
Due  to  the  hysteresis  of  the  JJs,  the  input  and  biasing  currents  must  both  be  reduced 
to  zero  to  allow  the  JJs  to  return  to  the  ZVS  and,  thus,  return  the  logic  gate  to  a 
logic  “low”. 

6.3.1  FD-TLM  implementation 

As  before,  the  circuit  is  placed  in  a  50x50x50  jum3  (x,  y ,  z)  box  filled  with  air  and  the 
walls  of  the  box  are  modeled  as  perfect  conductors.  The  bottom  surface  of  the  box 
is  used  as  the  superconducting  ground  plane  and  a  3  ym  thick  layer  of  SiC>2  is  placed 
at  the  bottom  of  the  box  on  which  the  rest  of  the  circuit  is  modeled.  All  conductors 
are  modeled  as  perfect  conductors  being  1  ym  thick,  and  generally  5  /an  wide.  Fig 
6.5.  shows  the  physical  layout  of  the  circuit. 


Figure  6.5:  Physical  layout  of  the  MVTL  circuit  used  in  FD-TLM  simulation. 
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Figure  6.6:  Results  of  simulation  of  the  MVTL  logic  gate. 

The  voltage  source  waveforms  are  shown  in  Fig.  6.6.  The  JJs  are  implemented 
using  the  same  parameters  as  before  with  the  exception  of  J2  having  a  critical  current 
and  junction  capacitance  three  times  that  of  Ji  and  J3. 

As  before,  the  JJs  are  implemented  as  a  node  centered  in  the  area  representing 
the  tunneling  barrier  of  the  JJ.  Appendix  D  lists  the  FD-TLM  data  set  describing 
the  MVTL  circuit.  The  simulation  is  performed  for  185  ps  with  a  1.7  fs  time  step 
dictated  by  the  1  fim  uniform  grid  spacing  used  in  the  FD-TLM  method.  FD-TLM 
method  results  are  shown  in  Fig.  6.6. 

6.3.2  Validation  of  the  FD-TLM  results 

To  perform  conventional  circuit  analysis  for  comparison  to  the  FD-TLM  results,  nodal 
analysis  is  performed  on  the  circuit  giving  the  following  equations: 
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Important  pieces  of  the  physical  layout  of  the  circuit  used  for  the  FD-TLM  method 
are  extracted  and  the  physical  dimensions  and  parameters  were  provided  to  the  Fas- 
tHenry  program  to  calculate  and  extract  the  inductances  and  mutual  inductances  as 
well  as  the  coupling  between  the  control  line  and  SQUID  loop  (see  Appendix  D). 
MAPLE  is  used  to  solve  for  a  set  of  first-order  differential  equations  (see  Appendix 
D)  and  these  equations  were  then  solved  using  the  fifth-order  time-adaptive  Runge- 
Kutta  method.  Results  of  this  method  are  shown  in  Fig.  6.6  along  with  the  FD-TLM 
method  results. 

Comparison  of  the  FD-TLM  method  results  and  the  Runge-Kutta  results  show  ex¬ 
cellent  agreement  therefore  validating  the  FD-TLM  method  for  modeling  the  MVTL 
circuit.  Furthermore,  the  results  for  both  the  FD-TLM  and  Runge-Kutta  methods 
agree  with  the  expected  logic  functioning  of  the  circuit  described  earlier.  Minor  differ¬ 
ences  can  be  attributed  to  inexact  values  for  the  inductances  and  magnetic  coupling 


calculated  with  FastHenry  as  will  as  parasitics  that  are  modeled  by  the  FD-TLM 
method  and  are  omitted  in  the  conventional  circuit  analysis. 


Chapter  7 
Conclusion 


With  the  advent  of  high  speed  digital  integrated  circuits,  analysis  and  design  is  be¬ 
coming  ever  more  difficult  with  the  need  for  determination  of  the  effects  of  propaga¬ 
tion  delay,  crosstalk,  dispersion,  signal  reflections,  signal  reflections  from  the  package 
walls,  radiation,  and  modal  dispersion.  Furthermore,  additional  parasitics  must  be 
analyzed  at  the  higher  operating  speeds  where  the  capacitances  and  inductances  often 
need  to  be  simulated  as  distributed  rather  than  lumped  components.  JJ  logic  circuits 
can  operate  at  frequencies  from  1  GHz  to  100  GHz  necessitating  consideration  of  high 
frequency  electromagnetic  effects  just  discussed. 

Quasi-static  simulation  methods  such  as  SPICE  are  inadequate  at  modeling  JJ 
logic  circuits  since  the  cross-sectional  dimensions  of  the  circuit  may  be  comparable 
to  the  wavelength  of  the  signals  involved.  Furthermore,  although  a  conventional 
circuit  simulation  method  can  analyze  propagation  delays,  crosstalk,  and  parasitics, 
such  analysis  requires  painstaking  effort  to  extract  the  values  of  capacitances  and 
inductances  within  the  circuit. 

The  FD-TLM  method  allows  an  individual  to  describe  the  physical  geometry  of  a 
J  J  logic  IC  where  the  the  parasitics  are  automatically  modeled  and  separate  extraction 
is  not  required.  The  FD-TLM  method  solves  Maxwell’s  curl  equations,  as  described 
in  Chapter  3,  allowing  simultaneous  simulation  of  the  operation  of  the  logic  gate 
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and  electromagnetic  parasitics.  As  a  result,  the  user  can  simulate  a  JJ  logic  circuit 
with  confidence  that  the  circuit  being  simulated  will  work  as  predicted.  The  only 
extraction  required  is  that  of  the  process  parameter  values  for  the  JJ  describing  the 
current- voltage  relationship  of  the  JJ. 

In  Chapters  5  and  6  several  different  logic  circuits  were  simulated  and  analyzed 
to  validate  the  FD-TLM  method  for  modeling  JJ  circuits.  Before  performing  conven¬ 
tional  circuit  analysis,  special  care  was  taken  in  designing  the  FD-TLM  layout  so  that 
parasitics  which  could  not  be  simulated  in  conventional  circuit  analysis  were  mini¬ 
mized.  This  often  required  several  trial  and  error  simulations  to  reduce,  for  example, 
the  lengths  of  conductors,  or  rearranging  the  geometry  of  the  layout. 

Upon  examining  the  FD-TLM  results,  earlier  statements  concerning  the  utility  of 
the  FD-TLM  method  were  substantiated.  That  is,  the  FD-TLM  method  often  pre¬ 
dicts  the  existence  of  parasitics  omitted  in  conventional  circuit  simulation.  Differing 
results  between  conventional  circuit  simulation  and  FD-TLM  simulation  usually  indi¬ 
cated  that  not  all  parasitics  had  been  included  in  the  conventional  circuit  simulation. 
Adjusting  the  conventional  circuit  simulation  to  take  these  effects  into  account  ulti¬ 
mately  lead  to  close  agreement  between  the  two  methods.  As  a  result,  the  FD-TLM 
method  often  became  a  tool  for  learning  details  about  the  JJ  and  JJ  logic  circuit 
operation  that  are  not  discussed  in  reference  papers  and  books. 

The  final  results  of  Chapters  5  and  6  showed  that  the  FD-TLM  method  is  not 
only  reliable  for  simulating  JJ  logic  circuits,  but  that  it  is  invaluable  as  well.  In  this 
thesis,  effort  was  focused  on  simulating  and  validating  the  results  of  simulation  for 
JJ  logic  circuits.  The  conclusion  was  that  with  proper  incorporation  of  both  dc  and 
ac  characteristics  of  the  JJ  device,  any  type  of  JJ  circuit  can  be  simulated  correctly. 

Simulations  of  JJ  logic  circuits  using  the  FD-TLM  method  took  approximately 
two  hours  on  average.  On  the  other  hand,  conventional  circuit  simulation  using  the 
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Runge-Kutta  method  often  took  approximately  15  minutes.  The  justification  for  using 
the  FD-TLM  method  in  this  situation  is  that  the  FD-TLM  method  provides  the  full 
electromagnetic  behavior  of  the  circuit  instead  of  using,  possibly  unjustified,  quasi¬ 
static  approximations.  In  addition,  the  process  required  to  perform  a  nodal  analysis 
for  each  circuit  and  extract  the  parasitic  inductance  and  capacitance  was  far  more 
time  consuming  than  the  two  hour  simulation  time  for  FD-TLM  simulation.  With  this 
new  method  for  simulating  Josephson  junction  circuits,  it  is  possible  for  fabrication 
procedures  to  be  directly  converted  to  an  FD-TLM  data  set  for  full  electromagnetic 
simulation  of  the  fabricated  JJ  IC  circuit.  This  means  that  the  FD-TLM  method  can 
be  used  to  directly  simulate  a  fabricated  IC  and  determine  problems  with  design  as 
well  as  theoretical  circuit  operation. 

Future  work  in  this  area  includes  extension  of  the  FD-TLM  method  to  model 
the  skin  effect  in  superconductivity  where  resistance  exists  at  high  frequencies  within 
superconductors.  The  skin  effect  becomes  important  at  frequencies  above  100  GHz, 
and  has  already  been  modeled  with  the  FD-TD  method  [18].  Future  work  may  also 
include  simulation  of  several  types  of  JJ  logic  circuits  appearing  in  print  as  well  as 
simulation  of  chip  packages.  Ultimately,  work  should  be  pursued  in  collaboration  with 
individuals  performing  fabrication  of  JJ  IC  circuits  to  compare  FD-TLM  simulation 
with  real-time  results  to  further  validate  the  method. 


Appendix  A 

FORTRAN  Code  for  JJ 
Implementation 


This  section  contains  the  portions  of  the  FORTRAN  source  code  for  the  FD-TLM 

method  which  was  updated  to  model  the  Josephson  junction  (JJ)  and  SQUID. 

C  PROGRAM  FDVCUFOR 

C  Modified  March,  1995  by  Christopher  G.  Sentelle 
C 

C  Josephson  Junctions  are  specified  by  the  following  parameters: 

C  JJIC,  the  critical  current,  JJPO,  the  plasma  oscillation  frequency 

C  JJG1,  subgap  conductance  (linear  resistance),  JJI1,  JJG2,  JJVS,  and  JJVT, 

C  parameters  based  on  manufacture  of  the  Josephson  Junction. 

C  In  order  to  include  magnetic  effects  on  the  critical  current, 

C  we  specify  PHIO. 

C  The  equation  describing  non-SC  conductance  is  taken  from 
C  Rollins,  Greg  J.,  "Numerical  Simulator  for  Superconducting 
C  Integrated  Circuits",  IEEE  Trans,  on  Comp. -Aided  Design, 

C  Vol  10,  No  2.,  p  246,  Feb.  1991. 

C 

C  The  Josephson  Junction  is  specified  by  ’R*  as  the  model  indicator 
C  and  indicates  implementation 

C  of  changing  conductance  specified  by  the  Josephson  Junction  model  at 

C  the  xyz  location  specified  by  the  >R)  statement. 

C 

C  In  addition  to  implementing  the  simulation  of  the  JJ, 

C  this  version  of  the  FDTLM  program  allows  coupling  of  the 

C  SQUID  device  by  forcing  a  couple  between  two  JJ’s  that  are 

C  connected  in  a  superconductive  ring.  This  method  alters  the 

C  critical  current  of  the  JJ’s  based  on  the  magnetic  flux  evident 

C  in  the  loop.  The  following  device  card  is  used  to  specify  a 

C  Josephson  Junction  SQUID. 

C 

C  H  XI  X2  YPLN  Z1  Z2  JJ1  JJ2 
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c 

C  where  X1,X2,Z1,  and  Z2  specify  the  area  of  the  superconducting  ring. 

C  JJ1  and  JJ2  are  the  number  specifiers  for  the  JJ's  that  are  affected 

C  by  the  superconductive  ring  and  are  coupled.  YPLN  specifies  the 
C  location  of  the  flux  area 
C 
C 

C  H  2  JJ  SQUID  coupling 

C  R  Josephson  Junction 

C 
C 

C  Read  DATA  for  the  Josephson  Junction  Model 
C  Josephson  Junctions  (R  or  r) . 

C 

C 

C  Get  the  X,Y,Z  directions  for  the  Josephson  Junction 
C 

1000  CALL  GETDAT(IXLO,FJUNK, A JUNK, ICURSR, AINPUT, 1) 

IF  ( (IXLO  .LT.  1)  .OR.  (IXLO  .GT.  NX))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8,482)  IXLO, NX 
STOP 
END  IF 
C 

CALL  GETD AT (IXHI , F JUNK , A JUNK , ICURSR , AINPUT , 1 ) 

IF  ((IXHI  .LT.  1)  .OR.  (IXHI  .GT.  NX))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8 , 482)  IXHI,  NX 
STOP 
END  IF 
C 

CALL  GETD AT (IYP0S , F JUNK , A JUNK , ICURSR , AINPUT , 1 ) 

IF  ( (IYP0S  .LT.  1)  .OR.  (IYP0S  .GT.  NY))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8 ,483)  IYP0S,NY 
STOP 
END  IF 
C 

CALL  GETD AT ( I ZL0 , F JUNK , A JUNK , I CURSR , AINPUT , 1 ) 

IF  ( (IZL0  .LT.  1)  .OR.  (IZLO  .GT.  NZ))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE (8, 484)  IZLO,  NZ 
STOP 
END  IF 
C 

CALL  GETD AT ( I ZHI , F JUNK , A JUNK , I CURSR , AINPUT , 1 ) 

IF  ((IZHI  .LT.  1)  .OR.  (IZHI  .GT.  NZ))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8 , 484)  IZHI,  NZ 
STOP 
ENDIF 
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c 

C  Get  the  Plasma  Frequency  Po 
C 

CALL  GETD AT ( I JUNK , J JPO , A JUNK , I CURSR , AINPUT , 2 ) 

C 

C  Get  the  Critical  Current  of  the  Josephson  Junction 
C 

CALL  GETDAT(I JUNK , JJIC , A JUNK , I CURSR , AINPUT ,2) 

C 

C  Get  the  sub-gap  conductance  for  the  Josephson  Junction 
C 

CALL  GETD AT ( I JUNK , J JG 1 , A JUNK , I CURSR , AINPUT , 2 ) 

C 

C 

CALL  G  ETD  AT ( I JUNK , J  J 1 1 , A JUNK , I CURSR , A INPUT , 2 ) 

C 

CALL  GETD AT ( I JUNK , JJG2 , A JUNK , ICURSR , AINPUT , 2 ) 

C 

CALL  GETD AT ( I JUNK , JJVS , A JUNK , ICURSR , AINPUT ,2) 

C 

CALL  GETD AT ( I JUNK , J J VT , A JUNK , I CURSR , AINPUT , 2 ) 

C 

CALL  GETD AT ( I JUNK , J JC J , A JUNK , ICURSR , AINPUT , 2 ) 

C 

CALL  GETD AT ( I JUNK , JJDEP , A JUNK , ICURSR, AINPUT ,2) 

C 

NJJ  =  NJJ  +  1 
C 

C  Check  whether  the  array  will  overflow. 

C 

IF  (NJJ  .GT.  LJJ)  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8 , 930)  LJJ 

930  FORMAT ( ’  The  number  of  Josephson  Junctions  ’  , 

1  ’is  greater  than’ /’ the  array  size  LJJ  =  ’ , 

2  15, Increase  LJJ  and  recompile.  * ) 

STOP 

END  IF 
C 

IJJ(1 ,NJJ)  =  (IXHI+IXL0)/2 
I J J(2 ,NJ J)  =  IYPOS 
I J J (3 ,NJJ)  =  (IZHI+IZLO) /2 
IJJ(4,NJJ)  =  IXLO 
I JJ (5 ,NJJ)  =  IXHI 
I JJ (6 ,NJJ)  =  IZLO 
IJJ(7,NJJ)  =  IZHI 

FJJ ( 1 ,NJ J)  =  JJPO 
FJJ(2 ,NJJ)  =  JJIC 
FJJ (3 , NJJ)  =  JJG1 
FJJ (4 ,NJ J)  =  JJI1 


C 
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FJJ(5,NJJ)  =  JJG2 
FJJ(6,NJJ)  =  JJVS 
FJJ(7 ,NJJ)  =  JJVT 
FJJ (8,NJJ)  =  JJCJ 
FJJ(IO.NJJ)  =  JJDEP 
FJJ ( 13 ,NJJ )  =  ZERO 
9999  F0RMAT(G10 .4) 

9998  F0RMAT(I5) 

C 

C 

GOTO  1500 
C 
C 

C  Data  Entry  for  the  Squid  Coupling  Device.  Parameter  ’H’ 
C 


1200  NSqUID  =  NsquiD  +  1 

IF  (  NSqUID  .GT.  LSqUID  )  THEN 
WRITE(8, 160)  AINPUT 
WRITE(8, 1201)  LSqUID 

1201  FORMAT(’  The  number  of  Squids  exceeds  the  maximum  allowed  ’,15) 

STOP 

END  IF 

CALL  GETDAT ( IXLO , FJUNK , AJUNK , ICURSR , AINPUT , 1 ) 

CALL  GETD AT ( I XHI , FJUNK , AJUNK , I CURSR .AINPUT , 1 ) 

CALL  GETDAT(IYPOS , FJUNK , AJUNK , ICURSR, AINPUT , 1) 

CALL  GETDAT(IZLO, FJUNK, AJUNK, ICURSR, AINPUT, 1) 

CALL  GETDAT(IZHI , FJUNK , AJUNK , ICURSR .AINPUT , 1 ) 

CALL  GETDAT ( J J 1 , FJUNK , AJUNK , I CURSR , AINPUT , 1 ) 

CALL  GETDAT ( JJ2 , FJUNK , AJUNK , ICURSR , AINPUT , 1 ) 

C 

C  Check  to  be  sure  there  are  no  obvious  errors  in  the  the  device  parameters 

C 

IF  ((IXLO  .LT.  1)  .OR.  (IXLO  .GT.  NX))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8,482)  NX 
STOP 
END  IF 


C 


C 


IF  ((IXHI  .LT.  1)  .OR.  (IXHI  .GT.  NX))  THEN 
WRITE (8, 160)  AINPUT 
WRITE(8 ,482)  NX 
STOP 
END  IF 


IF  ((IYPOS  .LT.  1)  .OR.  (IYPOS  .GT.  NY))  THEN 
WRITE(8, 160)  AINPUT 
WRITE(8 ,482)  NY 
STOP 
END  IF 


C 

IF  ((IZLO  .LT.  1)  .OR.  (IZHI  .GT.  NZ))  THEN 
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WRITE(8 , 160)  AINPUT 
WRITE(8 ,484)  NZ 
STOP 
END  IF 
C 

IF  ((IZHI  .LT.  1)  .OR.  (IZHI  .GT.  NZ))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8 ,484)  NZ 
STOP 
END  IF 
C 

IF  ((JJ1  .GT.  NJJ )  .OR.  (JJ2  .GT.  NJJ))  THEN 
WRITE(8 , 160)  AINPUT 
WRITE(8 , 1202) 

1202  FORMAT( J Specify  JJs  before  specifying  SQUID  couplings ) 
STOP 
END  IF 
C 

C  Place  the  information  into  the  Squid  array 
C 

ISQUID(1,NSQUID)  =  IXLO 
ISQUID (2 , NSQUID)  =  IXHI 
I SQUID (3 , NSQUID)  =  IYPOS 
ISQUID (4 , NSQUID )  =  IZLO 
I SQUID (5 , NSQUID)  =  IZHI 
ISQUID(6, NSQUID)  =  JJ1 
ISQUID(7, NSQUID)  =  JJ2 
FJJ ( 13 , J Jl)  =  ZERO 
FJJ (13 , J J2)  =  NSQUID 
C 

GOTO  1500 
C 
C 

C  *****  E  field  calculations 
C 

C$DOACROSS  LOCAL(I , J ,K) 

DO  5350  K  =  1,  NZ 
DO  5300  J  =  1,  NY 
DO  5250  I  =  1,  NX 

EX(I , J ,K)  =  EXE(I , J ,K)*EX(I , J ,K)  +  EXH(I,J,K)* 

1  (HZ(I , J ,K)  -  HZ(I , J-l ,K)  +  HY ( I , J , K- 1 )  -  HY(I,J,K)) 

C 

EY(I , J ,K)  =  EYE ( I , J , K ) *EY ( I , J , K )  +  EYH(I,J,K)* 

1  (HX(I,J,K)  -  HX(I,J,K-1)  +  HZ(I-1,J,K)  -  HZ(I,J,K)) 

C 

EZ(I , J ,K)  =  EZE(I,J,K)*EZ(I,J,K)  +  EZH(I,J,K)+ 

1  (HY(I , J ,K)  -  HY(I-1,J,K)  +  HX(I , J-l ,K)  -  HX(I,J,K)) 

5250  CONTINUE 

5300  CONTINUE 
5350  CONTINUE 
C 


c 

c 


Special  FD-TLM  Ey  field  equations  are  used  for  the  JJ. 


DO  5400  M=l,  NJJ 
I  =  IJJ(1,M) 

J  =  IJJ(2,M) 

K  =  IJJ(3,M) 

C 

C  Calculate  the  old  Ey  electric  field  from  the  new  Ey. 

EY(I , J ,K)  =  (EY(I, J,K)  -  EYH(I,J,K)* 

1  (HX(I , J ,K)-HX(I , J,K-1)+HZ(I-1 , J,K)-HZ(I , J ,K) )) 

2  /EYE(I, J,K) 

C 

C  Calculate  the  JJ  phase 

IF  (FJJ(13,M)  .EQ.  ZERO)  THEN 

FJJ(ll.M)  =  FJJ(11,M)  -  FJJ(1,M)*DELTAT*EY(I,J,K) 

ELSE 

IXLO  =  ISqUID(l,FJJ(l3,M)) 

IXHI  =  ISqUID(2,FJJ(13,M) ) 

IYPOS  =  ISQUID(3,FJJ(13,M) ) 

IZLO  =  ISqUID(4,FJJ(13,M) ) 

IZHI  =  ISqUID(5,FJJ(13,M) ) 

JJ1  =  ISqUID(6,FJJ(13,M) ) 

JJ2  =  ISqUID(7,FJJ(13,M)) 

FLUXJJ  =  ZERO 
DO  5500  J1  =  IXLO,  IXHI 
DO  5501  K1  =  IZLO,  IZHI 

FLUXJJ  =  FLUX J J+U0*HY ( J 1 , IYPOS ,K1 ) *DELTAL*U( Jl) 

1  *W(K1)/V( IYPOS) 

5501  CONTINUE 

5500  CONTINUE 

FLUXJJ  =  TWO*PI*FLUXJJ/FLUXO 
FJJ(11,M)  =  FJJ(11, JJ1)  -  FLUXJJ 
END  IF 

C  Prevent  phase  overflow. 

IF  (FJJ ( 11 ,M)  .GT.  (TWO+PI) )  FJJ(11,M)  =  FJJ(11,M)  -  TWO*PI 
C 

C  Calculate  the  total  flux  magnitude  at  the  JJ  node 

FLUXJJX  =  ZERO 

FLUXJJZ  =  ZERO 

C  DO  5401  I1=IJJ(4,M) ,IJJ(5,M) 

C  FLUXJJZ  =  FLUXJJZ  +  HZ(I1,J,K) 

C  5401  CONTINUE 

C  DO  5402  II  =  IJJ(6,M) ,IJJ(7,M) 

C  FLUXJJX  =  FLUXJJX  +  HX(I,J,Ii) 

C  5402  CONTINUE 

C  FLUXJJZ  =  FLUXJJZ  *  UO  *  FJJ(10,M) 

C  FLUXJJX  =  FLUXJJX  *  UO  *  FJJ(10,M) 

C  FLUXJJ  =  SqRT(FLUXJJZ*FLUXJJZ+FLUXJJX*FLUXJJX) 

C  IF  (FLUXJJ  .LT.  1.0D-200)  THEN 

C  FLUXJJ  =  ONE 

C  ELSE 
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C  FLUXJJ  =  ABS(SIN(PI*FLUXJJ/FLUXO)/ ( (PI*FLUXJJ)/FLUXO) ) 

C  END  IF 

FLUXJJ  =  ONE 
C 
C 

C  Calculate  the  JJ  superconducting  current  including 
C  critical  current  effects  from  the  magnetic  flux 
FJJ(12,M)  =  FJJ(2,M)*SIN(FJJ(11,M) )*FLUXJJ 
C 

C  Calculate  the  new  Ey  field  including  the  effect  of  the  JJ. 

EY(I,J,K)  =  EY(I,J,K)  *  (EYE (I , J ,K) 

1  -  HALF  *  EYH(I,J,K)  *  FJJ(2,M) 

2  *  C0S(FJJ(11 ,M) )  *  FJJ ( 1 ,M) 

3  *  DELTAT) 

4  +  EYH(I,J,K)  *  (HX(I,J,K)  -  HX(I,J,K-l) 

5  +  HZ(I-1,J,K)  -  HZ(I,J,K) 

6  +  FJJ(12 ,M) ) 

C 

5400  CONTINUE 
C 

C  Update  the  Josephson  Junction  conductances,  capacitances  remain  constant 
C 

C  Remember  the  following  assignments 

C  IJJ(1,I)  =  x  position  of  JJ 

C  IJJ(2,I)  =  y  position  of  JJ 

C  I JJ (3,1)  =  z  position  of  JJ 

C 

C  FJJ (1,1)  =  Plasma  Frequency  Po 

C  FJJ(2,I)  =  Critical  Current  IC 

C  FJJ(3,I)  =  sub-gap  conductance  G1 

C  FJJ(4,I)  =11  (parameter  based  on  manufacture) 

C  FJJ (5, I)  =  G2  (parameter  based  on  manufacture) 

C  FJJ (6, I)  =  VS  (parameter  based  on  manufacture) 

C  FJJ (7, I)  =  VT  (parameter  based  on  manufacture) 

C  FJJ (8, I)  =  Cj  Capacitance  of  the  JJ 

C  FJJ (9, I)  =  node  admittance  without  the  Josephson  Junction 

C  FJJ (10, I)  =  Depth  of  Josephson  Junction 
C  FJJ(11 ,1)  =  phi  for  the  JJ 
C  FJJ(12, I)  =  IJ  for  the  JJ 

C  FJJ(13,I)  =  Boolean,  Update  phase  based  on  voltage?  Used  for  Squid  couple 
C 

DO  6000  I  =  1 ,NJJ 
C 

C  Calculate  the  voltage  across  the  JJ. 

VJJ  =  -EY(IJJ(1, I) , IJJ(2,I) , IJJ(3,I) ) 

C 

GJR  =  FJJ(3,I)*VJJ+(FJJ(4,I)+FJJ(5,I)*ABS(VJJ))*(l/(l+ 

1  EXP ( (F J J ( 6 , I ) -VJJ ) /FJJ (7,1)))- 

2  1/ ( 1+EXP ((FJJ(6,I)+VJJ)/FJJ(7,I)))) 

C 


C  Calculate  conductance  to  give  required  current  at  specified  V 
IF  (VJJ  .EQ.  ZERO)  THEN 
GJJ  “  ZERO 
ELSE 

GJJ  =  ZO+GJR/VJJ 
END  IF 
C 

YSJ  =  TW0*H*FJJ(8,I)/(E0*DELTAL) 

EYE(IJJ(1,I),IJJ(2,I),IJJ(3,I))  =  ONE  - 
1  TWO*GJJ/ (GJ J+FJ J(9 , I)+YSJ) 

EYH(IJJ(1 , I) , IJJ(2,I) , IJJ(3, I) )  =  TWO  *  ZO 
1  /(GJJ+FJJ(9,I)+YSJ) 

C 

6000  CONTINUE 
C 

C  Calculations  for  coupling  between  JJ  devices  to  form  a  two  junction  SQUID 
C 

C  ISQUID(1 , I)  =  IXLO 

C  ISQUID(2,I)  =  IXHI 

C  ISQUID(3,I)  =  IZLO 

C  ISQUID(4, I)  =  IZHI 

C  ISqUID(5,I)  =  1st  Josephson  Device  to  be  linked 

C  ISQUID(6,I)  =  2nd  Josephson  Device  to  be  linked 

C 

C  DO  6001  I  =  1 ,NSQUID 

C  IXLO  =  ISQUID(1,I) 

C  IXHI  =  ISQUID(2,I) 

C  IYPOS  =  ISQUID(3,I) 

C  IZLO  =  ISQUID(4,I) 

C  IZHI  =  ISQUID (5 , I) 

C  JJ1  =  ISQUID (6 , I) 

C  JJ2  =  ISQUID(7 , I) 

C  FLUXJJ  =  ZERO 

C  We  first  calculate  the  flux,  Hy  going  through  the  area 
C  DO  6002  J  =  IXLO,  IXHI 

C  DO  6003  K  =  IZLO,  IZHI 

C  FLUXJJ  =  FLUXJ J+U0*HY( J , IYPOS , K) *DELTAL*U( J) *W (K) /V ( IYPOS) 

C  6003  CONTINUE 
C  6002  CONTINUE 

C  We  will  link  the  squid  JJ’s  with  the  pi*f lux/f luxo  relationship 

C  by  adjusting  the  phase  of  the  second  JJ  listed,  JJ2 

C  In  a  coupled  squid,  the  only  way  to  change  the  value  of 

C  the  phase  of  the  second  JJ  is  by  altering  the  first  JJ  or 

C  altering  the  flux.  This  altercation  is  only  performed  at  this 

C  stage.  In  the  main  iterative  loop,  the  phase  for  JJ2  of  this  couple 

C  is  not  allowed  to  change.  This  routine  will  set  its  phase. 

C 

C  The  next  step  will  be  to  determine  the  phase  differences 
C 

C  FLUXJJ  =  TWO*PI*FLUXJJ/FLUXO 

C 


C  We  now  adjust  the  coupled  Josephson  Junctions 
C  FJJ(11 , JJ2)  =  FJJ ( 11 , JJ1)+FLUXJJ 

C 

C  6001  CONTINUE 


Appendix  B 

JAWS  Simulation  Data 

The  section  contains  the  FD-TLM  Data  Set,  the  FastHenry  extraction  data  set, 
MAPLE  V  results,  and  the  C  program  for  performing  the  conventional  circuit  simu¬ 
lation  for  the  JAWS  circuit. 


B.l  FD-TLM  Data  Set 

NFRCJL5 

T  RCJL  Josephson  Junction  Logic 

♦Generated  by  FDTGRAPH,  copyright  1993,  by  Christopher  G.  Sentelle 
♦Data  in  format  to  be  used  by  FD-TLM  copyright  by  Dr.  Robert  H.  Voelker 
♦University  of  Nebraska-Lincoln 
♦ 

♦Modified  for  simtime  and  pulses  on  Sept  6,  1995 
♦ 

♦Created  on:  Thu  Jun  29  13:02:33  1995 
♦ 

♦  In  this  simulation,  we  try  to  reduce  parasitic  capacitance  and  inductance 

♦  with  a  couple  of  methods.  First,  all  of  the  conducting  lines  are 

♦  made  to  be  infinitely  thin  to  take  care  of  parasitic  inductance. 

♦  There  should  be  little  to  no  crosstalk  in  this  circuit  because 

♦  everything  is  orthogonalized  as  much  as  possible.  We  also  reduce 

♦  the  area  of  the  JJs  in  order  to  try  to  reduce  a  capacitance  that  may 

♦  be  occuring  between  the  junction  overriding  the  capacitance  we  are 

♦  trying  to  create  via  the  Josephson  Junction  itself.  We  are  also 

♦  moving  the  entire  circuit  at  a  higher  level,  2  microns  from 

♦  ground  in  order  to  reduce  the  parasitic  capacitance  to  ground. 

♦ 

♦Medium  material  used  throughout 
♦Relative  permittivity 
E  1  50  1  25  1  50  1  1  1 
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♦Conductivity  throughout 
L  1  50  1  25  1  50  0  0  0 
♦Magnetic  susceptibility 
M  1  50  1  25  1  50  0  0  0 
♦Relative  permeability 
U  1  50  1  25  1  50  1  1  1 
♦Add  for  a  Si02  substrate 
E  1  50  121  50  3. 53. 53. 5 

♦  Josephson  Junction  #1 

R  6  7  3  12  13  3 . 039e+15  7e-05  0.004  0  0.1  0.002  0.0001  5e-13  5e-ll 

*  Josephson  Junction  #2 

R  6  7  3  22  23  3.039e+15  0.0001  0.004  0  0.1  0.002  0.0001  5e-13  5e-ll 
♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  8  3  3  3  8  -1  -2  -1 
L  993338-2-2-1 
L  683399-1-2-2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  12  14  3  3  3  8  -1  -2  -1 

L  15  15  3  3  3  8  -2  -2  -1 

L  12  14  3  3  9  9  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  9  11  3  3  6  8  -1  -2  -1 
L  12  12  3  3  6  8  -2  -2  -1 
L  9  11  3  3  9  9  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  253368-1-2-1 
L  663368-2-2-1 
L  253399-1-2-2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  8  3  3  9  11  -1  -2  -1 

L  9  9  3  3  9  11  -2  -2  -1 

L  6  8  3  3  12  12  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  7  3  3  12  13  -1  -2  -1 

L  8  8  3  3  12  13  -2  -2  -1 

L  6  7  3  3  14  14  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  18  20  4  4  3  18  -1  -2  -1 

L  21  21  4  4  3  18  -2  -2  -1 

L  18  20  4  4  19  19  -1  -2  -2 
♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  17  4  4  16  18  -1  -2  -1 
L  18  18  4  4  16  18  -2  -2  -1 
L  6  17  4  4  19  19  -1  -2  -2 
♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  8  4  4  14  15  -1  -2  -1 

L  9  9  4  4  14  15  -2  -2  -1 

L  6  8  4  4  16  16  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  7  4  4  12  13  -1  -2  -1 

L  8  8  4  4  12  13  -2  -2  -1 
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L  6  7  4  4  14  14  -I  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  2  5  4  4  16  18  -1  -2  -1 

L  6  6  4  4  16  18  -2  ~2  -1 

L  2  5  4  4  19  19  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  8  4  4  19  21  -1  -2  -1 

L  9  9  4  4  19  21  -2  -2  -1 

L  6  8  4  4  22  22  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  6  7  3  4  22  23  -1  -2  -1 

L  8  8  3  4  22  23  -2  -2  -1 

L  6  7  3  4  24  24  -1  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Conductor 
L  1  5  3  3  22  24  -1  -2  -1 

L  6  6  3  3  22  24  -2  -2  -1 

L  1  5  3  3  25  25  -1  -2  -2 

♦  These  values  of  resistance  had  to  be  corrected  by  hand. 
♦Infinitely  Thin  object,  xz-plane  Resistor  45  ohm  z-directed 
L  12  14  3  3  2  2  -2  -2  5555.56 

L  15  15  3  3  2  2  -2  -2  5555.56 

♦Infinitely  Thin  object,  xz-plane  Resistor  45  ohm  z-directed 
L  683322-2-2  5555.56 

L  993322-2-2  5555.56 

♦Infinitely  Thin  object,  xz-plane  Resistor  1  ohm 
L  113368  250000  -2  -2 
L  113399  250000  -2  -2 

♦Infinitely  Thin  object,  xz-plane  Resistor  12k  ohms 
L  18  20  4  4  2  2  -2  -2  20.8333 

L  21  21  4  4  2  2  -2  -2  20.8333 

♦Infinitely  Thin  object,  xz-plane  Resistor  45  ohm  x-directed 
L  1  1  4  4  16  18  5555.56  -2  -2 

L  1  1  4  4  19  19  5555.56  -2  -2 

♦  1  ps  =  600  iterations 

♦ 

♦  Minimum  Grid  Spacing 
A  le-06 

* 

♦Simulation  time  (ps)  300 
S  210000 

♦Backup  Interval  (ps)  200 
B  220000 

♦Plot  Interval  (ps)  0.1 
P  60 
* 

♦  A  Voltage  Source 

*  A  Pulse  waveform 

♦  Zero  Initial  Time  0(ps) 

♦  Rise  Time  15(ps) 

*  On  Time  185(ps) 

*  Fall  Time  15(ps) 


*  On  Voltage  -0.7 

*  Off  Voltage  0 

* 

VP  18  21  441  1Z0  9000  111000  9000  -0.7  0 

*  A  Voltage  Source 

*  A  Pulse  waveform 

*  Zero  Initial  Time  35(ps) 

*  Rise  Time  15(ps) 

*  On  Time  30(ps) 

*  Fall  Time  15(ps) 

*  On  Voltage  -0.002 

*  Off  Voltage  0 

* 

VP693311Z  21000  9000  57000  9000  -0.002  0 
V  P  12  15  3  3  1  1  Z  60000  9000  18000  9000  -0.002  0 
♦Voltage  Paths 
W  Y  1  3  19  4 
W  Y  1  2  13  4 
W  Y  1  2  7  4 
W  Y  1  2  7  10 
W  Y  1  3  7  17 
W  Y  -1  3  7  21 
* 

♦Current  Loops 
J  Z  17  21  3  5  4 
J  Z  11  15  2  4  4 
J  Z  5  9  2  4  4 

J  Z  5  9  2  4  10 

J  X  3  5  15  19  4 

J  X  2  4  21  25  3 

*  Variable  Mesh  Array 

*  Variable  Mesh  in  the  X  direction 
G  X  1  50  1 

*  Variable  Mesh  in  the  Y  direction 
G  Y  1  25  1 

*  Variable  Mesh  in  the  Z  direction 
G  Z  1  50  1 
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B.2  FastHenry  Data  Set 


This  is  the  data  set  provided  to  the  FastHenry  program  to  determine  the  values  of 
the  two  parasitic  inductances  in  the  JAWS  circuit. 

*  Determination  of  inductances  in  the  JAWS  circuit  for 

*  verification  of  simulation  at  home 

*  Performed  on  18  July  1995 

*  Uses  fasthenry 
.Units  urn 

*  Use  a  high  value  for  conductivity,  Superconductor 
.Default  nhinc  =  1  nwinc  =  5  sigma=1.0e20  z=2.5  w=3  h=0.1 

*  Ground  Plane 

gl  xl=0  yl=0  zl=0 
+  x2=50  y2=0  z2=0 

+  x3=50  y3=50  z3=0 

+  segl=20  seg2=20 
+  thick=l 

*  Setup  for  the  system 
N1  x=3  y=7 

N2  x=24  y=7 
N3  x=22 . 5  y=6 
N4  x=22 . 5  y=2 
♦Connect  the  segments 
El  N1  N2 
E2  N3  N4 

♦Make  needed  electrical  connections  between  segments. 

.equiv  N2  N3 
♦Define  output 
.external  N1  N4 

.freq  fmin=ie9  fmax=le9  ndec=l 
.  end 
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B.3  MAPLE  V  Results 

These  are  the  results  after  using  MAPLE  V  to  solve  the  set  of  differential  equations 
obtained  from  nodal  analysis  in  terms  of  a  set  of  first  order  derivatives  of  each  variable. 
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>  eqnl:=(Vl(t)-Va(t) )/Rin+(Vl(t)-Vb(t) )/Rin+Cd*diff(Vl(t) ,t)+Vl(t) 

>  /Rd+Il ( t )  =  0:  1  ' 

eqnI:=mizim ,  VHO-WO  lC/a, . Who . 

Rin  Rin 


,( d  \ 


+  +I1<,)=0 


>  eqn2 : =-11 ( t ) +Cj  *di f  f ( V2 ( t ) -V3 ( t ) , t ) +G1 ( V2 ( t ) -V3( t) )+Iclwsin(phil 

>  ( t ) ) =0 ; 

eqn2  :=-Il(/)  +  Q  ^|v2(/)]-[|v3(/)jj  +  Gl(V2(/)-V3(0)  +  /c7  sin(phil(O)  =  0 

>  eqn3 : =(V3 ( t )-Vso( t ) )/Rso+Cj *di f f ( V3 ( t ) . t ) +G1 ( V3 ( t ) ) +Ic2*sin ( phi2 

>  J^J:C^diff(V2(t)-V3(t).t)-Gl(V2(t)-V3(t) )  -Icl*sin  ( phil  ( t )  )  +12  ( 

C<}n3  :=  ~  ^Rso  S0<^  +  Cj  (  J/  V3^ 1 ) )  +  G 1  (  V3(  ‘)  > +  Ic2  sin(  phi2(  t ) ) 


-  Q [[a  V2(0 j-^V3(/) JJ-  Gl(  V2(0  -  V3(0) -  Id  sin(phil(/))  + 12(0  =  0 
>  eqn4:=-I2(t)+Cout*diff(V4(t) . t ) +V4 ( t )/Rout=0 ; 

eqn4  :=  -12(0  +  Cout  ( —  "A‘ •  -V4^  -  « 


(|v+o) 


+  — ^-=0 
Rout 


>  eqn5 : =di£f (phil ( t ) . t ) =Po* ( V2 ( t ) -V3 ( t ) ) ; 


eqn 5  :=  Ft ph  ‘ 1  ( 0  =  Po  ( V2(  0  ~  V3(  0  ) 

>  eqn6 : =di f f ( phi2 ( t ) . t ) =Po*V3 ( t ) ; 

eqn 6  :=  ^  phi2(  0  =  T’o  V3(/) 

>  eqn7:=Vl( t)-V2( t)  =  Ll'di f f ( II ( t ) , t ) ; 

eqn?  :=Vl(t)-V2(t)  =  LI  (|ll(o) 

>  eqn8 : =V3 ( t ) -V4 ( t )  =  L2*di f f ( 12 ( t ) . t ) ; 


eqn8  :=  V3(0  -  V4(/)  =  L2  12(0  J 

>  ans : =solve( {eqnl . eqn2 . eqn3 . eqn4 , eqn5 , eqn6 . eqn7 . eqn8 } , (diff (VI (t) 

>  . tj.diff (V2(t) .t) .diff (V3(t) ,t) ,diff(V4(t) ,t) .diff (Il(t) ,t) .diff 

>  ( 12 ( t ). t) .diff ( phil ( t) . t) .diff (phi2(t) . t) }) : 

am  :=  [ft phl  1(/)  =  Po  ( V2(/)  ~  V3(0),  Jt  phi2(0  =  Po  V3  (/), 

ly./fV  2  RdVl(t)-  Rd'Va(t)-  RJVb(t)  +  V\(t)  Rin  +  U(t)  Rin  Rd  d . 

*  {)  CdRinRd - >jV2(t)  =  -( 

-2  11(0 /bo  +  V3(0~  Vso(0  +  Gl(V3(0)  Rso  +  Ic2  sin(phi2(0)  Rso  +  \2{t)  Rso 
+  Gl(  V2(0  -  V3(0)  Rso  +  IcI  sin(phil(0)  Rso) /(Cj Rso), £  V3 (/)  = 

-1 1  ( Q  Rso  +  V3(  Q  -  Vso(  Q  +  G 1  ( V3(  Q )  Rso  +  Ic2  sin(  phi2(  Q )  Rso  + 12(  Q  Rso 


Cj  Rso 
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3,,m  I2(t)  Rout  -V4(t)  a  -Vl(f)  +  V2(Q  a  V3(/)-V4(Q| 

aV4(?)=  coutR^  -a? 11(0  -  u  ’a/140'  12  [ 

>  ans[l]; 

|phil(/)=JPo(V2(0-V3(0) 

>  ans[2] ; 

|phi2(/)  =  POV3(0 

>  ans[ 3 ] ; 

a  2fo/Vl(/)-fo/Va(0-/^Vb(0  +  Vl(0/?>w  +  Il(Qi?fflfo/ 

a<V1(,)-'  CdRinRd 


> 


> 


> 


> 


ans[ 4]  ; 

|-V2(0  =  -(-2  I1(0£»>  +  V3(0-  Vso(/)  +  Gl(V3(0)^so  +  /c2sin(phi2(0)/iyo 

c/ 

+ 12(0  +  Gl(  V2(/)  -  V3(Q)  Rso  +  /c/  sin(phil(f))flso)/(Cyfoo) _ 

ans[5] ; 

£  -U(t)Rso  +  W3(t)- Vso( /)  +  G1(V3(/))  Rso  +  Ic2  sin( phi2( Q)foo  +  I2(Qfoo 

ftV3(Q— _ Cyfoo _ 

ans  [  6  ]  ; 


ans[7]  ; 


a  i£0£o£-V4(0 

a/  j  Cout  Rout 


-Vl(/)  + V2(/) 
LI 


V3Q)-V4(/) 

L2 


>  ans[8] : 
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B.4  C  Code 

This  is  a  portion  of  the  C  Code  used  to  solve  the  set  of  first  order  differential  equations 

provided  by  MAPLE  using  the  Runge-Kutta  fifth-order  method.  Conventional  circuit 

simulation  results  are  then  obtained. 

#include  "nrutil.c" 

#include  "odeint.c" 

#include  "rkqs.c" 

#include  "rkck.c" 

#include  "linint.c" 

#def ine  NOVAR  8 
#define  Icl  0.7e-4 
#define  Ic2  1 . Oe-4 
#def ine  Po  3.039ei5 
#def ine  Cd  0.292e-15 
#define  Cout  0.327e-15 
#def ine  LI  7.0e-12 
#def ine  L2  18.00e-12 
#def ine  C j 1  7.0e-13 
#define  Cj2  7.0e-13 
#define  Cj  5.0e-13 
#define  Rd  1.0 
#define  Rso  12.0e3 
#define  Rin  45.0 
#define  Rout  45.0 
#def ine  R1  45.0 
#def ine  PI  3.141592763 
#def ine  SIMTIM  350.0e-12 
#def ine  EPS  1.0e-2 


int  kmax=1200 ,kount ; 

float  *xp,**yp,dxsav=3. Oe-13; 


float  Gl (float  v) 

{ 

float  gl=4e-3; 
float  g2=0.1; 
float  il=0.0; 
float  vs=2e-3 ; 
float  vt=0.1e-3; 
float  ans; 


ans  =  gl*v  +  (il  +  g2*fabs(v) )*((l/(i+exp((vs-v)/vt) )) 
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-(l/(l+exp((vs+v)/vt) )) ) ; 


> 


return  ans ; 


float  Vso(float  t) 

{ 

float  starttime  =  0.0e-12; 
float  risetime  =  15.0e-12; 
float  ontime  =  185.0e-12; 
float  falltime  =  15.0e-12; 
float  hicur  =0.7; 


> 


if  (t  <=starttime) 
return  0.0; 

if  ((t>starttime)  &&  (t<=starttime+risetime) ) 
return  (t-starttime) *hicur/risetime ; 
if  ( (t>starttime+risetime)&&(t<=starttime+risetime+ontime) ) 
return  hicur; 

if ( (t>starttime+risetime+ontime)&&(t<=starttime 
+risetime+ontime+f allt ime) ) 

return  (hicur  -  (t-(starttime+risetime+ontime) ) 

*hi cur/falltime) ; 

else 

return  0.0; 


float  Va(float  t) 

{ 

float  starttime  =  35.0e-12; 
float  risetime  =  15.0e-12; 
float  ontime  =  95.0e-12; 
float  falltime  =  15.0e-12; 
float  hicur  =  2.0e-3; 

if  (t  <=starttime) 
return  0.0; 

if  ((t>starttime)  &&  (t<=starttime+risetime) ) 
return  (t-starttime) *hicur/riset ime; 
if  ( (t>starttime+risetime)&&(t<=starttime+risetime+ontime) ) 
return  hicur; 

if ( (t>starttime+risetime+ontime)&&(t<=startt ime+risetime+ontime 
+f alltime) ) 

return  (hicur  -  (t-(starttime+risetime+ontime) ) 

*hicur/f alltime) ; 

else 

return  0.0; 

} 

float  Vb(float  t) 

{ 

float  starttime  =  100.0e-i2; 
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float  risetime  =  15.0e-12; 
float  ontime  =  30.0e-12; 
float  falltime  =  15.0e-12; 
float  hicur  =  2.0e-3; 

if  (t  <=starttime) 
return  0.0; 

if  ((t>starttime)  &&  (t<=starttime+risetime) ) 
return  (t-startt ime ) *hicur/r iset ime ; 
if  ( (t>starttime+risetime)&&(t<=starttime+risetime+ontime) ) 
return  hicur; 

if ( (t>startt ime+riset ime+ont ime ) && ( t<=st artt ime+r iset ime+ont ime 
+f alltime) ) 

return  (hicur  -  (t-(starttime+risetime+ontime) ) 

*hicur/f alltime) ; 

else 

return  0.0; 


void  derive (float  t,  float  yv[] ,  float  dydt [] ) 
int  i ; 

/*In  this  new  system,  the  following  variables  apply: 
yv[l]=Vl(t) ; 
yv[2]=V2(t) ; 
yv [3] =phil(t) ; 
yv [3] =phi2(t) ; 

*/ 

float  Gl(float  v) ; 

float  Vso(float  t); 
float  Va(float  t); 
float  Vb(float  t); 


/*  The  following  set  of  equations  include  two  inductors,  one 
at  the  output  in  order  to  see  if  we  can  reduce  oscillation 
magnitude  in  the  output  when  the  circuit  is  in  HVS 

yvCl]  =  Vi(t) 
yv[2]  =  V2(t) 
yv  [3]  =  V3(t) 
yv  [4]  =  V4(t) 
yv [5]  =  Il(t) 
yv [6]  =  I2(t) 
yv[7]  =  phil(t) 
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yv[8]  =  phi2(t) 

*/ 

dydt[l]  =  ~(2.*Rd*yv[l]-Rd*Va(t)-Rd*Vb(t)+yv[l]*Rin 
+y v [5] *Rin*Rd)/(Cd*Rin*Rd) ; 

dydt[2]  =  -(-2. *yv[5] *Rso+yv[3]-Vso(t)+Gl(yv[3] )*Rso 
+Ic2*sin(yv[8] )*Rso 

+yv [6] *Rso+Gl (y v [2] -y v [3] ) *Rso+I cl*s in(yv [7] ) *Rso ) 

/ (Cj*Rso) ; 

dydt  [3]  =  “(-yv[5] *Rso+yv[3] -Vso (t)+Gl(yv [3]  )*Rso+Ic2*sin(yv [8] )*Rso 
+y v [6] *Rso)/ (Cj  *Rso) ; 

dydt  [4]  =  (yv[6]*Rout-yv[4] )/(Cout*Rout) ; 

dydt  [5]  =  -(-yv  [l]+yv  [2] )/Ll ; 

dydt  [6]  =  (yv  [3]  -y v  [4]  )  /L2 ; 

dydt  [7]  =  Po*(yv[2]-yv[3]  ) ; 

dydt  [8]  =  Po*yv[3]; 


int  main(){ 

FILE  *file; 
float  *vstart; 
float  result; 
int  i ; 

int  nok,  nbad; 
float  ts; 

/*  This  next  simulation  will  simulate  the  MVTL  circuit.  Simulation 
will  be  performed  with  this  system  followed  by  simulation  with  the 
FDTLM  method. 

The  following  configurations  will  be  used. 

We  shall  drive  a  current  through  the  system  that  is  just  under  the 
maximum  critical  current  of  2.0e-4  amps.  We  will  then 
apply  a  control  current  to  the  control  circuitry  linked 
to  our  JJ’s  through  the  inductance  L.  This  should  alter 
the  configuration  of  the  system  and  create  a  HVS. 


We  shall  analyze  the  input  current,  the  currents 
through  each  JJ,  and  the  voltage  across  each  JJ  to 
determine  what  state  the  system  is  in! 

Later  models  will  add  the  output  resistance  to  shunt  the 
current  along  with  an  input  current  through  a  voltage  and 
a  resistor.  An  inductance  will  be  used  later  to  model 
effects  of  the  loop  and  its  self  inductance. 

We  want  to  see  what  the  basic  operation  of  the  SQUID  should  be*/ 


printf ("Allocating  memory. An"); 
xp  =  vector(l ,kmax) ; 
yp  =  matrix (1 ,N0VAR, 1 ,kmax) ; 
vstart  =  vector(l ,N0VAR) ; 


/*Clear  the  Matrices*/ 
printf ("Clearing  memory\n"); 
f  or ( i= 1 ; i<=N0VAR ; i++ ) { 
vstart [i] =0.0; 

> 

printf ("Performing  Runge  Kutta  5th  order  adaptive  integration . \n" 
odeint (vstart , NOVAR ,0.0 , SIMTIM , EPS , 1 . Oe- 13 , HMIN , took , 
tobad , derive , rkqs ) ; 

printf  ("nok=  '/,d  nbad=  */fd\n"  ,nok,nbad) ; 

printf  ("Memory  allocated  =  7,fK\n" ,  (kmax*NOVAR*sizeof  (f  loat) 
+kmax*sizeof (float) )/1000.0) ; 

printf ("Simulation  over,  calculating  and  printing  results\n"); 


/♦Print  Va  results*/ 
f ile=f open("jrcjlva. dat" , "w") ; 
f or (ts=0 . 0 ; ts<=SIMTIM ; ts+=i . Oe-13) 

f  printf  (file,  n,/,g  */lg\n"  ,ts ,  Va(ts) ) ; 
f close(f ile) ; 

/♦Print  Vb  results*/ 
f ile=f open(" jrcjlvb .dat" , "w") ; 
f or (ts=0.0;ts<=SIMTIM;ts+=l. Oe-13) 

f  printf  (f  ile ,  "'/*g  %g\n"  ,ts ,  Vb(ts)  )  ; 
f close(f ile) ; 

/♦Print  output  results  VI*/ 
f ile=f open("jrcjlvl .dat" , "w") ; 
f or(i=l ; i<=kount ; i++) 

f  printf  (file ,  M'/,g\n"  ,yp  [1]  [i]  ) ; 


f close(f ile) ; 


/♦Print  out  the  output  current*/ 
printf(" Interpolating:  Output  current\nM); 
f ile=f open("jrcjlv2 .dat" , Mw") ; 
f or (ts=0 . 0 ; ts<=SIMTIM; t s+=l . 0e-13){ 

linint (&xp[l] ,&yp[4]  [1] ,kount ,ts ,&result) ; 
fprintf  (file,  M</.g  */,g\n"  ,ts ,  result) ; 

} 

f close(f ile) ; 


/♦Print  out  the  Vso*/ 
f ile=f open(" jrcjlvo .datH , "w") ; 
f or (ts=0 ; ts<=SIMTIM; ts+=l . Oe-13) 

fprintf  (file,  '7,g  ,/.g\n"  ,ts  ,Vso(ts)  )  ; 
f close(f ile) ; 

/♦Print  source  current*/ 
file=fopen(" jrcjlcin.dat"  ,  V)  ; 
for(i=l; i<=kount ; i++) 

fprintf  (file,"*/, g\n" ,  (Vso(xp  [i]  )-yp  [2]  [i]  )/Rso)  ; 
f close(f ile) ; 

/♦Print  the  time  scale  */ 

file  =  f open("timescal . dat" , "w") ; 

f or(i=l ; i<=kount ; i++) 

fprintf  (f  ile ,  M,/.g\nM  ,  xp  [i]  ) ; 
f close(f ile) ; 


/♦Determine  shunted  current*/ 
printf  ("Done . \n"); 


return(O) ; 


Appendix  C 

2  JJ  DC  SQUID  Simulation  Data 


The  section  contains  the  FD-TLM  Data  Set,  the  FastHenry  extraction  data  set, 
MAPLE  V  results,  and  the  C  program  for  performing  the  conventional  circuit  simu¬ 
lation  for  the  2  JJ  DC  SQUID  circuit. 

C.l  FD-TLM  Data  Set 

NDATATEST 

T  2  JJ  Squid  Simulation  8/7/95 

♦Generated  by  FDTGRAPH,  copyright  1993,  by  Christopher  G.  Sentelle 
♦Data  in  format  to  be  used  by  FD-TLM  copyright  by  Dr.  Robert  H.  Voelker 
♦University  of  Nebraska-Lincoln 
* 

♦ 

♦Created  on:  Mon  Aug  7  14:42:20  1995 
♦Medium  material  used  throughout 
♦Relative  permittivity 
E  1  50  1  50  1  50  1  1  1 

♦Conductivity  throughout 
L  1  50  1  50  1  50  0  0  0 
♦Magnetic  susceptibility 
M  1  50  1  50  1  50  0  0  0 
♦Relative  permeability 
U  1  50  1  50  1  50  1  1  1 

♦Add  a  Si02  Substrate  or  Ground  Plane  Insulator 
E  150  131  50  3. 53. 53. 5 

♦  Josephson  Junction  #1 

R  11  15  3  4  8  3 . 039e+15  0.0001  0.004  0  0.1  0.002  0.0001  5e-13  5e-ll 

♦  Josephson  Junction  #2 

R  11  15  3  25  29  3.039e+15  0.0001  0.004  0  0.1  0.002  0.0001  5e-13  5e-ll 

♦  Coupling  device 
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H  10  18  4  9  24  1  2 

*3  dimensional  object,  Superconductor 
L  5  9  4  4  4  29  -1  -1  -1 
*Top  Surface 

L  5  9  5  5  4  29  “1  -2  -1 

L  10  10  S  S  4  29  -2  -2  -1 

L  5  9  5  5  30  30  -1  -2  -2 

♦Right  Side  Surface 

L  5  9  4  4  30  30  -1  -1  -2 

L  10  10  4  4  30  30  -2  -1  -2 

♦Back  Side  Surface 

L  10  10  4  4  4  29  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  10  15  4  4  25  29  -1  -1  -1 

♦Top  Surface 

L  10  15  5  5  25  29  -1  -2  -1 

L  16  16  5  5  25  29  -2  -2  -1 

L  10  15  5  5  30  30  -1  -2  -2 

♦Right  Side  Surface 

L  10  15  4  4  30  30  -1  -1  -2 

L  16  16  4  4  30  30  -2  -1  -2 

♦Back  Side  Surface 

L  16  16  4  4  25  29  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  10  15  4  4  4  8  -1  -1  -1 

♦Top  Surface 

L  10  15  5  5  4  8  -1  -2  -1 

L  16  16  5  5  4  8  -2  -2  -1 

L  10  15  5  5  9  9  -1  -2  -2 

♦Right  Side  Surface 

L  10  15  4  4  9  9  -1  -1  -2 

L  16  16  4  4  9  9  -2  -1  -2 

♦Back  Side  Surface 

L  16  16  4  4  4  8  -2  -1  -I 

♦3  dimensional  object,  Superconductor 

L  17  22  4  4  25  29  -1  -1-1 

♦Top  Surface 

L  17  22  5  5  25  29  -1  -2  -1 

L  23  23  5  5  25  29  -2  -2  -1 

L  17  22  5  5  30  30  -1  -2  -2 

♦Right  Side  Surface 

L  17  22  4  4  30  30  -1  -1  -2 

L  23  23  4  4  30  30  -2  -1  -2 

♦Back  Side  Surface 

L  23  23  4  4  25  29  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  17  22  4  4  4  8  -1  -1  -1 

♦Top  Surface 

L  17  22  5  5  4  8  -1  -2  -1 

L  23  23  5  5  4  8  -2  -2  -1 

L  17  22  5  5  9  9  -1  -2  -2 


♦Right  Side  Surface 
L  17  22  4  4  9  9  -1  -1  -2 

L  23  23  4  4  9  9  -2  -1  -2 

♦Back  Side  Surface 
L  23  23  4  4  4  8  -2  -1  -1 
♦3  dimensional  object,  Superconductor 
L  19  22  4  4  9  24  -1  -1  -1 
♦Top  Surface 

L  19  22  5  5  9  24  -1  -2  -1 

L  23  23  5  5  9  24  -2  -2  -1 

L  19  22  5  5  25  25  -1  -2  -2 

♦Right  Side  Surface 

L  19  22  4  4  25  25  -1  -1  -2 

L  23  23  4  4  25  25  -2  -1  -2 

♦Back  Side  Surface 

L  23  23  4  4  9  24  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  11  18  2  2  4  8  -1  -1  -1 

♦Top  Surface 

L  11  18  3  3  4  8  -1  -2  -1 

L  19  19  3  3  4  8  -2  -2  -1 

L  11  18  3  3  9  9  -1  -2  -2 

♦Right  Side  Surface 

L  11  18  2  2  9  9  -1  -1  -2 

L  19  19  2  2  9  9  -2  -1  -2 

♦Back  Side  Surface 

L  19  19  2  2  4  8  -2  “1  “1 

♦3  dimensional  object.  Superconductor 

L  11  18  2  2  25  29  -1  -1  -1 

♦Top  Surface 

L  11  18  3  3  25  29  -1  -2  -1 

L  19  19  3  3  25  29  -2  -2  -1 

L  11  18  3  3  30  30  -1  -2  -2 

♦Right  Side  Surface 

L  11  18  2  2  30  30  -1  -1  -2 

L  19  19  2  2  30  30  -2  -1-2 

♦Back  Side  Surface 

L  19  19  2  2  25  29  -2  -1  -1 

*3  dimensional  object,  Superconductor 

L  17  18  3  3  4  8  -1  -1  -1 

♦Top  Surface 

L  17  18  4  4  4  8  -1  -2  -1 

L  19  19  4  4  4  8  -2  -2  -1 

L  17  18  4  4  9  9  -1  -2  -2 

♦Right  Side  Surface 

L  17  18  3  3  9  9  -1  -1  -2 

L  19  19  3  3  9  9  -2  -1  -2 

♦Back  Side  Surface 

L  19  19  3  3  4  8  -2  -1  -1 

*3  dimensional  object,  Superconductor 

L  17  18  3  3  25  29  -1  -1  -1 


♦Top  Surface 

L  17  18  4  4  25  29  -1  -2  -1 

L  19  19  4  4  28  29  -2  -2  -1 

L  17  18  4  4  30  30  -1  -2  -2 

♦Right  Side  Surface 

L  17  18  3  3  30  30  -1  -1  -2 

L  19  19  3  3  30  30  -2  -1  -2 

♦Back  Side  Surface 

L  19  19  3  3  28  29  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  1  4  4  4  18  17  -1  -1  -1 

♦Top  Surface 

L  1  4  6  5  15  17  -1  -2  -1 

L  5  5  5  5  15  17  -2  -2  -1 

L  1  4  5  5  18  18  -1  -2  -2 

♦Right  Side  Surface 

L  1  4  4  4  18  18  -1  -1  -2 

L  5  5  4  4  18  18  -2  -1  -2 

♦Back  Side  Surface 

L  5  5  4  4  15  17  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  23  34  4  4  11  22  -1  -1  -1 

♦Top  Surface 

L  23  34  5  5  11  22  -1  -2  -1 

L  35  35  5  5  11  22  -2  -2  -1 

L  23  34  5  5  23  23  -1  -2  -2 

♦Right  Side  Surface 

L  23  34  4  4  23  23  -1  -1  -2 

L  35  35  4  4  23  23  -2  -1  -2 

♦Back  Side  Surface 

L  35  35  4  4  11  22  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  30  34  4  4  3  10  -1  “1  -1 

♦Top  Surface 

L  30  34  5  5  3  10  -1  -2  -1 

L  35  35  5  5  3  10  -2  -2  -1 

L  30  34  5  5  11  11  -1  -2  -2 

♦Right  Side  Surface 

L  30  34  4  4  11  11  -1  -1  -2 

L  35  35  4  4  11  11  -2  -1  -2 

♦Back  Side  Surface 

L  35  35  4  4  3  10  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  35  49  4  4  14  22  -1  -1  -1 

♦Top  Surface 

L  35  49  5  5  14  22  -1  -2  -1 

L  50  50  5  5  14  22  -2  -2  -1 

L  35  49  5  5  23  23  -1  -2  -2 

♦Right  Side  Surface 
L  35  49  4  4  23  23  -1  -1  -2 

L  50  50  4  4  23  23  -2  -1  -2 


♦Back  Side  Surface 
L  50  50  4  4  14  22  -2  -1  -1 
*3  dimensional  object,  Superconductor 
L  19  22  6  6  3  39  -1  -1  -1 
♦Top  Surface 

L  19  22  7  7  3  39  -I  -2  -1 

L  23  23  7  7  3  39  -2  -2  “1 

L  19  22  7  7  40  40  -1  -2  -2 

♦Right  Side  Surface 

L  19  22  6  6  40  40  -1  -1  -2 

L  23  23  6  6  40  40  -2  “1  -2 

♦Back  Side  Surface 

L  23  23  6  6  3  39  “2  “1  -1 

♦3  dimensional  object,  Superconductor 

L  2  18  6  6  35  39  -1  -1  -1 

♦Top  Surface 

L  2  18  7  7  35  39  -1  -2  -1 
L  19  19  7  7  35  39  Back  Side  Surface 
L  35  35  4  4  11  22  -2  -1  -1 
♦3  dimensional  object,  Superconductor 
L  30  34  4  4  3  10  -1  -1  -1 
♦Top  Surface 

L  30  34  5  5  3  10  -1  -2  -1 

L  35  35  5  5  3  10  -2  -2  -1 

L  30  34  5  5  11  11  -1  -2  -2 

♦Right  Side  Surface 

L  30  34  4  4  11  11  -1  -1  -2 

L  35  35  4  4  11  11  -2  -1  -2 

♦Back  Side  Surface 

L  35  35  4  4  3  10  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  35  49  4  4  14  22  -1  “1  -1 

♦Top  Surface 

L  35  49  5  5  14  22  -1  -2  -1 

L  50  50  5  5  14  22  -2  -2  -1 

L  35  49  5  5  23  23  -1  -2  -2 

♦Right  Side  Surface 

L  35  49  4  4  23  23  -1  -1  -2 

L  50  50  4  4  23  23  -2  -1  -2 

♦Back  Side  Surface 

L  50  50  4  4  14  22  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  19  22  6  6  3  39  -1  -1  -1 

♦Top  Surface 

L  19  22  7  7  3  39  -1  -2  -1 

L  23  23  7  7  3  39  -2  -2  -1 

L  19  22  7  7  40  40  -1  -2  -2 

♦Right  Side  Surface 
L  19  22  6  6  40  40  -1  -1  -2 

L  23  23  6  6  40  40  -2  -1  -2 

♦Back  Side  Surface 


L  23  23  6  6  3  39  -2  -1  -1 
*3  dimensional  object,  Superconductor 
L  2  18  6  6  35  39  -1  “1  -1 
♦Top  Surface 

L  2  18  7  7  35  39  “1  -2  -1 
L  19  19  7  7  35  39  -2  -2  -1 
L  2  18  7  7  40  40  -1  -2  -2 
♦Right  Side  Surface 
L  2  18  6  6  40  40  -1  -1  -2 
L  19  19  6  6  40  40  -2  -1  -2 
♦Back  Side  Surface 
L  19  19  6  6  35  39  -2  -1  -1 

*  The  following  conductivity  values  were  modified  to 

*  inaccuracy  of  the  Fdtgraph  program. . 

♦Infinitely  Thin  object,  xz-plane  SRC  RES  60  OHMS 
L  30  34  4  4  2  2  -2  -2  2777.78 

L  35  35  4  4  2  2  -2  -2  2777.78 

♦Infinitely  Thin  object,  xz-plane  LOAD  RES  16  OHM 
L  50  50  4  4  14  22  6250  -2  -2 

L  50  50  4  4  23  23  6250  -2  -2 

♦Infinitely  Thin  object,  xz-plane  CNTR  SRC  RES  16  OHMS 
L  19  22  6  6  2  2  -2  -2  12500 

L  23  23  6  6  2  2  -2  -2  12500 

♦Infinitely  Thin  object,  xz-plane  CNTR  LOAD  RES  1  OHM 
L  1  1  6  6  35  39  166666.67  -2  -2 

L  1  1  6  6  40  40  166666.67  -2  -2 

*  1  ps  =  600  iterations 

* 

*  Minimum  Grid  Spacing 
A  le-06 

+ 

♦Simulation  time  (ps)  120 
S  111000 

♦Backup  Interval  (ps)  130 
B  120000 

♦Plot  Interval  (ps)  0.1 
P  60 

* 

♦  A  Voltage  Source 

♦  A  Pulse  waveform 

♦  Zero  Initial  Time  0(ps) 

♦  Rise  Time  5(ps) 

+  On  Time  100(ps) 

♦  Fall  Time  5(ps) 

♦  On  Voltage  -0.01 

♦  Off  Voltage  0 

♦ 

VP  30  35  441  1Z0  3000  60000  3000  -0.01  0 

♦  A  Voltage  Source 

*  A  Pulse  waveform 

+  Zero  Initial  Time  40(ps) 
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*  Rise  Time  5(ps) 

*  On  Time  20(ps) 

*  Fall  Time  5(ps) 

*  On  Voltage  -0.004 

*  Off  Voltage  0 

* 

VP  19  23  661  1Z  24000  3000  12000  3000  - 
♦Voltage  Paths 
W  Y  1  3  32  5 
W  Y  1  5  21  4 
W  Y  1  5  9  37 
W  Y  1  3  18  6 
W  Y  1  3  18  27 
W  Y  1  3  12  6 
W  Y  1  3  12  27 
W  Y  1  3  26  16 
W  Y  1  3  44  18 
* 

♦Current  Loops 
J  Z  29  35  3  5  5 
J  Z  18  23  5  7  4 
J  X  5  7  34  40  9 
J  X  3  8  3  9  18 
J  X  3  5  24  30  18 
J  X  3  5  3  9  12 
J  X  3  5  24  30  12 
J  X  3  5  10  23  26 
J  X  3  5  13  23  44 

♦  Variable  Mesh  Array 

♦  Variable  Mesh  in  the  X  direction 
G  X  1  50  1 

♦  Variable  Mesh  in  the  Y  direction 
G  Y  1  50  1 

♦  Variable  Mesh  in  the  Z  direction 
G  Z  1  50  1 


0.004  0 


i 
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C.2  FastHenry  Data  Set 


This  is  the  data  set  provided  to  the  FastHenry  program  to  determine  the  values  of 
the  two  parasitic  inductances  in  the  2  JJ  DC  SQUID  circuit. 

♦Determination  of  Inductance  for  Squid  Loop 

*  August  7,  1995  ,  Christopher  Sentelle 

♦  Uses  fasthenry 
.Units  urn 

*  Use  a  high  value  for  conductivity,  Superconductor 
.Default  nhinc  =  1  nwinc  =  3  sigma=1.0e20  z=3  w=5  h=l 

*  Ground  Plane 

gl  xi=0  y 1=0  zl=0 
+  x2=55  y2=0  z2=0 

+  x3=S5  y3=55  z3=0 

+  segl=25  seg2=25 
+  thick=l 

*  Setup  for  the  system 

♦  Squid  Loop 

N5  x=29 . 0  y=3 . 0 
N6  x=31 . 0  y=5 . 5 
N7  x=40 . 0  y=5 . 5 
N8  x=42 . 5  y=3 . 0 
N9  x=42 . 5  y=15 . 0 
N10  x=42 . 5  y=17 . 0 
Nil  x=42 . 5  y=29 . 0 
N12  x=40 . 0  y=26 . 5 
N13  x=31 . 0  y=26 . 5 
N14  x=29 . 0  y=29 . 0 
N15  x=29 . 0  y=17 . 0 
N16  x=29 . 0  y=15 . 0 

*  Control  Line 

N17  x=29 . 0  y=2 . 0  z=5.0 
N18  x=29 . 0  y=39 . 0  z=5.0 
N19  x=31 . 0  y=36 . 5  z=5.0 
N20  x=48 . 0  y=36 .5  z=5.0 
♦Connect  the  segments 

♦  Squid  Loop 
E3  N16  N5  w=4 
E4  N6  N7 

E5  N8  N9 
E6  N10  Nil 
E7  N12  N13 
E8  N14  N15  w=4 

♦  Control  Line 
E9  N17  N18  w=4 
E10  N19  N20 

♦  Added  for  Mf  calculations 
Ell  N9  N10 


♦Make  needed  electrical  connections  between  segments. 

.equiv  N5  N6 

.equiv  N7  N8 

.equiv  Nil  N12 

.equiv  N13  N14 

.equiv  N18  N19 

♦Define  output 

♦Changed  for  Mf  calculation 

.external  N16  N15 

.external  N17  N20 

.freq  fmin=le9  fmax=le9  ndec=l 

.  end 


C.3  MAPLE  V  Results 


These  are  the  results  after  using  MAPLE  V  to  solve  the  set  of  differential  equations 
obtained  from  nodal  analysis  in  terms  of  a  set  of  first  order  derivatives  of  each  variable. 


>  eqnl:=Cj*dif f (Vl(t) . t)+Ic*sin(phil(t) )+G(Vl( t ) ) -II ( t ) =0 ; 

eqnl  :=  Cj  ^  Vl(oj  +  /c  sin(phil(0)  +  G(  Vl(f ))  - 11(0  =  0 

>  eqn2:=Cj’dif f (V2(t ) ,t)+Ic’rsin(phil(t)-2*PI*(Mf*I3(t)+L2*I2(t)-Ll 

>  *Il(t) )/FLo)+G(V2(t) )-I2(t)=0: 

eqn2  := 

o(|v2(0]w„in(phil(0-2n(^13(')t^(,)-t'l'('))),G(V2(,))-l2(,) 
=  0 


> 


> 


> 


eqn3 : =11 ( t ) +12 ( t ) -14 ( t ) +15 ( t ) =0 : 

eqn3  :=  11(0  + 12(0  -  14(0  + 15(/)  =  0 
eqn4 : =V4 ( t )/Rcout-I3 ( t ) =0 ; 


eqn4:=^i-  13(0  =  0 
H  Rcout  v 


eqn5 : =13 ( t )+(V5 ( t )-Vcin( t } )/Rcin=0 ; 

V5(Q- Vcin(Q 
Rein 


eqnS  I3(  t )  +  - 


=  0 


>  eqn6:=I4(t)+(V6(t)-Vso(t))/Rin=0: 

eqn6  :=  I4(  / )  + - - - =  0 

Rin 

>  eqn7:=-I5(t)+V7(t)/Rout=0; 

— <'>♦!?- 

>  eqn8:=V3(t)-Vl(t)=Ll*diff ( I 1 ( t ) . t ) -Ml*dif £ ( 13 ( t ) . t ) ; 

e9nS:=V3(0-V1(/)  =  L/^|ll(o]-A//j^|l3(o) 

>  eqn9 : =V3 (t)-V2(t)=L2*diff(I2(t) , t) +M2*di ff(I3(t).t); 

eqn9  :=  V3(  0  -  V2(  / )  =  1.2  (jj;  I2(  / )  j  +  M2  I3(  / )  j 

>  eqnlO : =V5 (t)-V4(t)=L3*diff(I3(t) . t ) -Ml *di f f ( II ( t ) , t ) +M2*di f f ( 12 (. 

eqnlO  :=  V5(/)  -  V4(/)  =  L3  (jt  I3(o] -  Ml  11(f)]  +  M2  I2(o] 

>  eqnll : =V6 ( t ) -V3 (t)=L4*diff(I4( t) . t) ; 

eqnl  1  :=  V6(0  -  V3(/)  =  L4  ^^14(0^ 

>  eqnl 2 : =V3 ( t ) -V7 (t)=L5*diff(I5(t) . t) : 

eqnl 2  :=  V3(  t )  -  V7(  0  =  L5  {j(  I5(  t )  j 

>  eqnl 3 : =di f f ( phil ( t ) , t)=Po*Vl(t) ; 

eqnl3  :=jtphi\(t)  =  l>o  V1(0 


> 

> 


> 


> 


> 


> 


> 


> 


> 


> 


> 


> 


> 


> 


> 


eqnlOa : =subs( {V5 ( t ) =solve ( eqn5 , V5 ( t ) ) , V4 ( t )=solve( eqn4 , V4 ( t ) ) } , e 
qnlO )  ; 

eqnlOa 

{ 13(f)  -  -~^]  Rein  -  13 (/)  Rcout  =  L3  I3(f)]-M  11(f)]  +  M2{ft  BO)] 

eqnlla : =subs( V6 ( t )=solve( eqn6 . V6 ( t ) ) , eqnll ) ; 

eqnlla  :=  1 14(f)  -  Rin  -  V3 (/)  =  L4  14(f)] 

eqnl2a : =subs ( V7 ( t ) =solve ( eqn7 , V7 ( t ) ) , eqnl2 ) ; 

eqnI2a  :=  V3(f )  -  15(0  Rout  =  L5 fj]  15(f)] 


eqnl: 

Q{iv,(0)+/csin(phi,(0)+G(V1(,))"I1(/)=0 

eqn2 ; 

Q(|v2(0),/c»{phiK0-2ni"m,H^2(0-"l'(,))).0(V2(0)-l2(0  = 
0 

eqn3 : 

_ I1(Q +  12(0 -14(0 +  15(0  =  0 _ 

eqn3 : 

V3(f )- Vl(f)  =  £/(|  11(f)]- A// 13(f)] 


eqn9; 

V3(f)-V2(f)  =  ^|l2(f)]  +  M2^|l3(f)] 

eqnlOa; 

{I3(0_:^)^rt-I3(0/?COM,=£5(lI3(,))~M/(ln(0)+iW2(lI2(0) 

eqnlla : 

eqnl2a ; 

V3(0  -  I5(  t)Rout  =  L5  I5(/)J 

eqn8a : =subs( V3 ( t ) =solve ( eqnl2a , V3 ( t ) ) , eqn8 ) ; 

eqn8a  :=  15(f)  Rout  +  U 15(f)]  -  Vl(/)  =  LI  (|ll(o)  -  Ml  I3(o] 

eqn9a : =subs( V3 ( t ) =solve ( eqnl 2a . V3 ( t ) ) , eqn9 ) ; 

eqn9a  :=  15(f)  Rout  +  LS  15(f)]  -  V2(f )  =  L2  [|  12(f)]  +  M2  13(f)] 

eqnl lb ; =subs( V3 ( t ) =solve ( eqnl 2a , V3 ( t ) ) , eqnlla ) ; 


eqnllb  ~  |l4(f)  -  ^°-]  Rin  - 15(f)  Rout  -  L5  I5(/)J  =  Z.V  14(f)] 

>  eqnl ;  ~~~  — 

Cy{ivl(0]+/csin(phil(0)+G(V1(0)' 11(0=0 

>  eqn2;  ~  ^  ~  ~ 

Q(gv2(,)),/c4hi1(,,-2nW^l)^/«>>-t'll<'»),0(V2(0)-l2(,) 

0 

>  e q n  3 ;  . . .  ~ . .  '  “  " 

11(f)  +  12(0  -I4(/)  +  15(0  =  0 

>  eqn8a;  -  _____  —  ~  ' 

15(0  Rout  +  L5  15(f)]  -  V1(0  =  U  (f  11(f)]  -  Ml  13(f)] 

>  eqn9a;  :  ~ 

I5(  f )  Rout  +  L5  15(f)]  -  V2(  f )  =  12  I2(  f )]  +  M2  [|  13(f)] 

>  eqnlOa;  ~~~  ~~  ~  ~ 

>  eqnllb;  ’ 

{ I4(  0  -  ^°-]  Rin  -15(f)  Rout  -  L5  I5(  t)j  =  L4  14(f)] 

>  eqnllc : =subs ( 14 ( t ) =solve ( eqn3 . 14 ( t ) ) .eqnllb) : 

eqnllc  := -^ll(f)  +  12(f)  +  15(f)  -  V*°(0]  Rin  - 15( f )  Rout  -  L5  ^  I5( f )]  = 

L4{j(  11(f)  +  12(f)  +  15(f)] 

>  eqnl; 

Q  (|  VI(f )]  +  Ic  sin(phil(f))  +  G(  Vl(f ))  -  H(f)  =  0 

>  eqn2;  ~ ~  ~ 

o(|v2(o)+/csi(phn(„  ^nattmumL IzMM »),0(v2(<>)  _  mo 
0 

>  eqn8a;  ~ 

I5(f)/?OWf  +  iLi^|l5(f)]-Vl(f)  =  L/^|u(f)]-M^|l3(f)] 

>  e  q  n  9  a ;  ”  . 

15(f)  Rout  +  L5  (f  15(f)]  -  V2(  f )  =  L2  I2(  f )]  +  M2  I3(  f )] 


>  eqnlOa; 


> 


> 


> 

> 

> 


{I3(/) -^Sr]  Rcin  ~n{l)  Rcout  = 13  (| 13^))  -  M}  (I  IK'))  +  M2  (|  I2(o) 

eqnllc;  ~  — — 

{ll(0  +  1 2(t)  + 15(0  -  - 15(0  Rout  -  L5  15(0  j  = 

eqnl3; 


^phii(0=PoVi(0 


ans : =solve( (eqnl . eqn2 . eqn8a . eqn9a, eqnlOa , eqnllc . eqnl3 } , (diff (VI ( 
t) .t) ,dif£(V2(t) .tj.diff (Il(t).t) .diff ( 12 ( t ) .t) ,diff(I3(t) . t) ,di 
ff(I5(t),t) ,diff (phil ( t ) . t ) } ) : 


|phii(0=^vi(0,|v2(0= 


Ic +  2  n nz,2  12(0-2  nz,/ 11(0  j _ G(  V2(0)  + 12(0 


Cj 


f  1 1  ( t )  =  (  Vso(  t)  M22  L5  +  L5  L4  M2  I3(  t )  Rcout -L4V2(t)L3  L5  -  L5  L4  M2  Vcin(  / ) 


+  L5  L4 L3  V1(0  +  L5  L4  M2 13(0  Rein  +  Ml  L3  L4 13(0  Rein  +  Ml  L5 L4 13(0  Rcout 

-  Ml  L5  M2  11(0  Rin  +  Ml  L5  M2  Vso(/)-  Ml  L5  M2 15 (t)Rin-  Ml  L5  M2  \2(t)Rln 

-  Ml  L5  L4  Vcin(0  -  L5  M22  VI (/)  +  I5(/)  Rin  L3  L5  L2  +  L5  L3  VI (/)  L2 

-II  (t)Rin  M22  L5  +  11(0  Rin  L3  L5  L2  - 12(0  Rin  M22  15  + 15  13(0  Rcout  Ml  L2 
+  L5 13(  t )  Rcin  Ml  L2  - 15 M2  Ml  V2(0  -  L5  Vcin(/)  Ml  L2  -  Vso(0  L3  L5  L2 

+ 12(0  Rin  L3  L5  L2  -  15(0  Rin  M22  L5  +  MI  L4  M2 15 (t)Rout  -  M22  L4  Vl(/) 

-  M2  Ml  L4  V2  ( 0  +  M22  L4  I5(  t)  Rout  +  L4 13(t)  Rcin  Ml  L2  -  L4  Vcin(  t)  Ml  L2 

+  L4  I3( 0  Rcout  Ml  L2  -  L4  L3 15(0  Rout  L2  +  L4  L3  Vl(/)  L2)/( %  1 ),  ^  I3( 0  =  -  ( 

-L5  L4 13(0  Rcout  L2  -  L5  L4 13(0  Rcin  L2  +  M2  LI  L5  V2(t)  -  M2  LI  Vso(t)  L5 
+  L5  L4  Vcin(0  L2  +  M2  LI  12(0  Rin  L5  +  M2  LI  15(0  Rin  L5  +  M2  LI  11(0  Rin  L5 
-M2LIL4 15( 0  Rout  +  M2  LI  L4  V2( /)  +  L4  Ml  15(0  Rout  L2  +  L5  Ml  Vso( t)  L2 
+  L5  M2  L4  V2(0~  L5  Ml  II (/)  Rin  L2  -  L5  Ml  15(0  Rin  L2  +  LI  L5  L4  Vcin(/) 

-  V1(0  Ml  L5  L2  -  V1(0  Ml  L4  L2  -  Vl(/)  MI  L5  L4  -  LI  L5  L4 13(0  Rcin 

-  LI  L5 L4 13(0  Rcout  +  L1  L4  Vein (/)  L2  -  LI  L5 13(0  Rcout  L2  -  LI  L5 13(0  Rcin  L2 
+  LI  L5  Vein  (t)L2  -  LI  L4 13(0  Rcout  L2  -  LI  L4 13(0  Rcin  L2-L4M2  VI  (/)  L5 

+  L4  Ml  V2(0  L5  -  L5  Ml  12(0  Rin  L2)/{% 1 ),  1 12(/)  =  (m!2  L5  Vso (0 

-  Ml2  V2(0  L5  -  L5  L4  M2  13(0  Rcout  +  L5 L3  LI  U(t)  Rin  +  L4  V2(  t)L3L5 


-  LS  L3  LI  Vso(/)  +  L5  L4  M2  Vein (0  +  V2 (/)  L3  LI  L5  +  M2  LI  L5  Vein (/) 

-  L5  L4L3V\{t)  +  L5  L3  LI  12 (/)  Rin  -  L5  L4  M2 13(0  Rein  +  L5  L3 LI  I5(/)  Rin 

-  M2  LI  L5 13(0  Reout  -  M2  LI  L5 13(0  Rein  -  Ml 2  L5 11(0  Rin  -  Ml 2  L5 12(0  Rin 

-  Ml2  L5 15(0  Rin  -  Ml  L5L4 13(0  Rein  -  Ml  LS  L4 13(0  Reout  -Ml  L5  M2 11(/)  Rin 
+  Ml  L5  M2  Vso(0  -  Ml  L5  M2 15(0  Rin  -  Ml  M2  Vl(/)  L5  -  Ml  LS  M2 12(0  Rin 

+  Ml  LS  L4  Vcin(0  -  Ml 2  L4  V2(t)  +  Ml  L4  M2  15(0  Rout  +  Ml2  L4  [5(0  Rout 

-  Ml  L4  M2  V1(0 -L3  LI  L4 15(0  Rout  +  L3  LI  L4  V2 (/)  -  L4 13(0  Rein  M2  LI 

+  L4  Vcin(0  M2  LI  -  L4 13(0  Reout  M2  Ll)/(%  1  ),J^I5(0  =  {m!2  L4  V2 (/) 

-  2  Ml  L4  M2 15(0  Rout  -  Ml2  L4 15(0  Rout  +  Ml2  Vso(0  L2  +M22  LI  Vso(t) 

+  M2 2  L4  V1(0  +  Ml  L4 M2  V\{t)  +  L3  LI  L4 15(0  Rout  -  L3  LI  L4  V2(0 

-  M22  LI  11(0  Rin  -  M22  LI  12(0  Rin  -  M22  LI  15(0  Rin  -  M22  LI  15(0  Rout 
+  M2  Ml  L4  V2(0  -  M22  L4 15(0  Rout  -  Ml2  11(0  Rin  12  ~  Ml 2  12(0  Rln  L2 

-  Ml2  15(0  Rin  L2-M1 2  15(0  Rout  L2  +  L3  LI  11(0  Rin  L2  +  L3  LI  12(0  Rln  12 
+  L3  LI  15(0  Rin  L2  -  L3  LI  Vso(  t)L2  +  L3  LI  15(  t)  Rout  L2  +  L4 13(  t )  Rein  M2  LI 

-  L4 13(0  Rein  Ml  L2  -  L4  Vcm(t  )M2  LI  +  L4  Vcm(t)Ml  L2  +  L4 13(0  Reout  M2  LI 

-  L4  13(0  Reout  Ml  L2  +  L4 1.3  15(0  Rout  L2  -  L4  L3  V1(0  L2)/(%  1 ), 

<LV\n\-  /csin(phil(0)  +  G(Vl(0)-Il(Q  , 
dt  ()  Cj 

%1  :=  Ml 2  LS  L4  +  2  Ml  LS  L4  M2  -  L3  LI  LS  L4  +  M22  LI  LS  +  M2 2  LI  L4+  M2 2  LS  L4 

+  Ml2  LS L2  +  Ml 2  L4  L2  -  L3  LI  LS  L2  -  L3  LI  L4  L2  -  L3  LS  L4  L2 

>  ans[l];  ~~  ”  7“  — 

|phii(0  =  ^vi(0 

>"ans[2] ;  ”  “ 

fv2(,)- 

f„i^i‘(')^^ri^l^n^i2(,)-2ni/n(,)j_0(v2(>))t[,(0 

Q 

>  ans[3]  ;  '  ~ — 

§  11(0  =  (  Vso( 0  M22  LS  +  LS  L4  M2  13(0  Reout  -  L4  V2(0  L3  LS  -  LS  L4  M2  Vein (/) 

+  LS  L4  L3  V1(0  +  LS  L4  M2 13 (t)Rcin  +  Ml  LS L4 13(0  Rein  +  Ml  LS  L4 13(0  Reout 
-Ml  LS  M2 11(0  Rin  +  Ml  LS  M2  Vso (/)  -  Ml  LS  M2  15(0  Rin  -  Ml  LS  M2 12(0  Rin 

-  Ml  LS  L4  Vcin(0  -  LS  M2 2  Vl (/)  + 15(0  Ri»L3  LS  L2  +  LS  L3  Vl {t)L2 


- 1 1(0  Rin  M22  L5  +  11(0  Rin  L3  L5  L2-Y2(t )  Rin  M22  L5  +  L5 13(0  Rcout  Ml  L2 
+ L5 13(f)  Rein  Ml  L2  -  L5  M2  Ml  V2(/)  -  L5  Vein  (t)Ml  L2  -  Vso(/)  L3  L5  L2 

+ 12(0  R‘n  13  L5  L2  - 15(0  Rin  M2 2  L5  +  Ml  L4  M2 15(0  Rout  -  M21  L4  Vl(/) 

-  M2  Ml  L4  V2(0  +  M2 2  L4 15(0  Rout  +  L4 13(0  Rein  MI  L2-L4  Vcin(f)  Ml  L2 
+  L4 13(0  Rcout  Ml  12  -  L4  L3 15(0  Rout  L2+L4 13  V1(0  L2)\(m12L5L4 

+  2M1 L5L4M2-L3L1L5L4  +  M22  LI  L5+M22L1 L4  +  M22L5L4  +  M12L5L2 

+  Ml 2  L4L2-L3  LI  L5  L2-L3LIL4L2-  L3  L5  L4  L2) 

>  an s  [  4  ]  ;  ~~~~~  — - 

J*I3(0  =  -  (-L5  L-f  13(0  Rcout  L2  -  L5  L4 13(0  RcinL2  +  M2L1 L5  V2 (/) 

-M2  LI  Vso (t)L5  +  L5 L4  Vcin(/)  L2  +  M2LlY2(t)  Rin  L5  +  M2L1 15(0  R'n ^ 

+  M2  LI  11(0  Rin  L5  -  M2  LI  L4 15(0  Rout  +  M2  LI  L4V2(t)+  L4  Ml  15(0  Rout  L2 
+  L5  Ml  Vso(0  L2  +  L5  M2  L4  V2(0  - 15 Ml  11(0  Rin  L2  -  L5  Ml  15(0  Rin  L2 
+  LI  L5  L4  Vcin(0  -  V1(0  Ml  L5  L2  -VI (/)  Ml  L4  L2  -  VI (/)  Ml  L5  L4 

-  LI  L5  L4 13(0  Rein  -  LI  L5  L4 13(0  Rcout  +  LI  L4  Vcin(0  L2  -  LI  L5 13(0  Rcout  L2 

-  LI  L5 13(0  Rein  L2  +  LI  L5  Vcin(/)  L2  -  LI  L4 13(0  Rcout  L2  -  LI  L4 13(0  Rein  L2 

-  L4  M2  V1(0  L5  +  L4  Ml  V2 (/)  L5  -  L5  Ml  12(0  Rin  L2)/{m1 2  L5L4  +  2  Ml  L5  L4  M2 

-  L3  LI  L5  L4  +  M22L1 L5  +  M22L1  L4  +  M22  L5  L4  +  M1 2  L5  L2  +  M1 2  L4  L2 
-L3  LI  L5  L2  -  L3  LI  L4  L2  -  L3  L5  L4  L2) 

>  ans[5] ;  —  ~ 

^  12(0  =  (.Ml2 L5  Vso(0  -Ml2  V2(0  L5 - L5  L4 M2  13(0  Rcout  +  L5  L3  LI  11(0  Rin 

+  L4  V2(0  L3  L5  -  L5  L3  LI  Vso(0  +  L5L4M2  Vein (/)  +  N2(t)L3  LI  L5 
+  M2L1L5  Vcin(0  -  L5  L4  L3  W\(t)  +  L5  L3  LI  Y2(t)Rin  -L5L4M2 13(0  Rein 

+  L5  L3  LI  15(0  Rin  -  M2  LI  L5 13(0  Rcout  -M2L1L5 13(0  Rein  -  Ml 2 1511(0  Rin 

-Ml2  L5 12(0  Rin  - Ml 2  L5 15(0  Rin  -  Ml  L5  L4 13(0  Rein  -  Ml  L5  L4  13(0  Rcout 

-  MI  L5 M2 11(0  Rin  +  Ml  L5 M2  Vso(0  -MIL5M2 15(0  Rin -Ml  M2  VI (/)  L5 

-M1L5M2 12(0  Rin  +  Ml  L5  L4  Vein (/)  -  Ml 2  L4  V2 (/)  +  Ml  L4  M2 15(0  Rout 
+  M12  L4 15(0  Rout  -  Ml  L4  M2  VI (/)  -  L3  LI  L4 15(0  Rout+L3  LI  L4  V2(0 

-  L4  13(0  Rein  M2  LI  +  L4  Vcin(0  M2L1-L4 13(0  Rcout  M2  Ll)/(M12  L5  L4 

+  2  MILS  L4  M2  -L3  LI  L5  L4  +  M22  LI  L5  +  M22  LI  L4  + M22  L5  L4  + Ml2  L5  L2 

+  M12  L4L2-  L3  LI  L5  L2-L3L1L4L2-  L3  L5  L4  L2) 

>  ans[6 ] :  ”  ' 

jt  I5(  0  =  (Ml 2  L4  V2( 0  -  2  Ml  L4  M2 15(  t )  Rout  -  Ml2  L4l5(t)  Rout  +  Ml2  Vso(  t )  L2 
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+  /W22  LI  \so{t)+M22  L4  Vl(/)  +  MI  L4  M2  VI (/)  +  L3  LI  L4 15(0  Rout 

-  L3  LI  L4  V2 (t)-M22 LI  U(t)Rin-M22  LI  12(0  Rin- M2 2  LI  15(0  Rin 

-  M22  LI  15(0  Rout  +  M2  MI  L4  V2 (t)-  M22  L4 15(0  Rout  -  MI 2 11(0  Rin  L2 

-  Ml 2 12( 0  Rin  L2  -  MI 2 15( 0  Rin  L2-  Ml2 15( t)  Rout L2  +  L3LI  U(t)RinL2 
+  L3  LI  12(0  Rin L2  +  L3  LI  I5(  t)  Rin  L2  -  L3 LI  Vso (/)  L2  +  L3  LI  15(0  Rout  L2 

+  L4 13(0  Rein  M2  LI  -  L4 13(t)  Rein  MI  L2  -  L4  Vcin(0  M2  LI  +  L4  Vcin(r)  Ml  L2 

+  L4 13(0  Rcout  M2L1-L4 13(0  Rcout Ml  L2  +  L4 L3 15(0  Rout  L2  -  L4 L3  VI (/)  Z.2) /( 
Ml2  L5  L4  +  2  Ml  L5  L4  M2  -  L3  LI  L5  L4  +  M22  LI  L5  +  M22  LI  L4  +  M22  L5  L4 
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C.4  C  Code 


This  is  a  portion  of  the  C  Code  used  to  solve  the  set  of  first  order  differential  equations 

provided  by  MAPLE  using  the  Runge-Kutta  fifth-order  method.  Conventional  circuit 

simulation  results  are  then  obtained. 

#include  "nrutil.c" 

#include  "odeint.c" 

#include  "rkqs . c" 

#include  "rkck.c" 


#def ine 

NOVAR  6 

#def ine 

L3 

6. Oe-12 

#def ine 

L4 

6. Oe-12 

#def ine 

Ml 

2. Oe-12 

#def ine 

M2 

2. Oe-12 

#def ine 

Mf 

6.068e-12 

#def ine 

LI 

8 . 39e-12 

#def ine 

L2 

8 . 39e-12 

#def ine 

Ic 

1.0e-6 

#def ine 

G1 

4 . Oe-3 

#def ine 

11 

0.0 

#def ine 

G2 

0.1 

#def ine 

Vs 

2.  Oe-3 

#def ine 

Vt 

1.0e-4 

#def ine 

Io 

1.0e-4 

#def ine 

Po 

3 . 039elS 

#def ine 

Cj 

5 . Oe-13 

#def ine 

Rein  16.0 

#def ine 

Rcout  1.0 

#def ine 

Rin  60.0 

#def ine 

Rout  16.0 

#def ine 

FLo  2 . 068e-15 

#def ine 

PLT  0.1e-12 

#def ine 

PI 

3.141592763 

#def ine 

SIMTIM  120 . Oe-12 

#def ine 

EPS  1.0e-2 

int  kmax=1300 ,kount ; 

float  *xp,**yp,dxsav=1.0e-13; 


float  G(float  v) 

{ 


float  ans; 


97 


ans  =  Gl*v  +  (II  +  G2*f abs (v) )*( (l/(l+exp( (Vs-v)/Vt ) ) ) 

- ( 1/ ( 1+exp ( ( Vs+v ) /Vt ) ) ) ) ; 

return  ans ; 


float  Vso(float  t) 

{ 

float  starttime  =  0.0e-12; 
float  risetime  =  5.0e-12; 
float  ontime  =  lOO.Oe-12; 
float  falltime  =  5.0e-12; 
float  hicur  =  0.010; 

if  (t  <=starttime) 
return  0.0; 

if  ( (t>starttime)  kk  (t<=starttime+risetime) ) 
return  (t-starttime) *hicur/risetime ; 
if  ( (t>starttime+risetime)&&(t<=starttime+risetime+ontime) ) 
return  hicur; 

if  (  ( t >  s t  art t  ime+r  i  s e t  ime+ont  ime )  kk  ( t  <=s t ar tt  ime+r  i  s  et  ime 
+ontime+f alltime) ) 

return  (hicur  -  (t-(starttime+risetime+ontime) ) 

*hicur/f alltime) ; 

else 

return  0.0; 

’  > 

float  Vcin(float  t) 

{ 

float  starttime  =  25.0e-12; 
float  risetime  =  5.0e-12; 
float  ontime  =  20.0e-12; 
float  falltime  =  5.0e-12; 
float  hicur  =  0.003; 

if  (t  <=starttime) 
return  0.0; 

if  ((t>starttime)  kk  (t<=starttime+risetime) ) 
return  (t-starttime) *hicur/ risetime ; 
if  ( (t>startt ime+risetime)&&(t<=starttime+risetime+ontime) ) 
return  hicur; 

if ( (t>startt ime+r i set ime+ont ime ) kk (t <=startt ime+r iset ime+ont ime 
+f alltime) ) 

return  (hicur  -  (t-(starttime+risetime+ontime) ) 

*hicur/f alltime) ; 

else 

return  0.0; 
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void  derive(float  t,  float  yv[] ,  float  dydt  [] ) 

{ 

int  i ; 

/*In  this  new  system,  the  following  variables  apply: 
y[l]=Vl 
y  [2]  =V2 
y  [31=11 
y  [4]=I2 
y  [5] =phil 

*/ 

float  G(float  v); 

float  Vso(float  t); 

float  Vcin(float  t); 


dydt [l]=-(Io*sin(yv[6] )+G(yv[l] )-yv[3] )/Cj ; 

dydt [2]=(Io*sin((-yv[6] *FLo+2*PI*(Mf *yv[5]+L2*yv[4] -Ll*yv [3] ))/FLo)- 
G(yv [2] )+yv[4] )/Cj ; 

dydt [3] =-(M2*M2*Rout*Vso(t)+M2*Ml*Rout*Vso (t )-M2*Ml*Rout*yv [3] 
*Rin-M2*Ml*Rout*yv [4] *Rin-M2*M2*yv [1] *Rin-M2*M2*yv [1] 
*Rout+L3*L2*yv [l] *Rin+L4*L2*Rout*yv[4] *Rin- 
L4*L2*Rout*Vso(t)+L4*L2*Rout*yv [3] *Rin+L4*L2*yv[l]*Rin 
+L4*L2*yv  [l] *Rout-L3*L2*Rout*Vso(t )+L3*L2*Rout*yv [3] 
*Rin+L3*L2*Rout *yv [4] *Rin+L3*L2*yv [1] *Rout-M2*M2*Rout 
*yv [3] *Rin-M2*M2*Rout*y v [4] *Rin-M2*Ml*yv [2] *Rin 
-M2*Ml*yv [2] *Rout-Vcin(t )*L2*Ml*Rin-Vcin(t)*L2*Ml*Rout 
+yv [5] *Rcout*L2*Ml*Rin+yv [5] *Rcout*L2*Ml*Rout+yv [5] *Rcin 
*L2*Ml*Rin+yv [5] *Rcin*L2*Ml*Rout)/ 

( (L2*L1*L4+L2*L1*L3-L2*M1*M1-L1*M2*M2 ) * (Rin+Rout ) ) ; 
dydt [4] =-(-Ml*Ml*yv[2]*Rin+Ml*Ml*Rout*Vso(t)-Ml*Ml*yv[2] *Rout 

+Ll*L4*Rout*yv[3] *Rin-Ll*L4*Rout*Vso(t)+Ll*L4*Rout*yv[4] 
*Rin+Ll*L4*yv  [2] *Rin+Ll*L4*yv [2] *Rout 
-Ll*L3*Rout*Vso(t)+Ll*L3*Rout*yv [3] *Rin+Ll*L3*Rout 
*yv  [4]  *Rin+Ll*L3*yv  [2]  *Rin+Ll*L3*yv  [2]  *Rout-Ml*Ml*Rout 
*yv [3] *Rin-Ml*Ml*Rout*yv [4] *Rin+M2*Ml*Rout*Vso (t ) 
-M2*Ml*Rout*yv  [3] *Rin-M2*Ml*Rout*yv[4] *Rin-M2*Ml*yv [1] *Rin 
-M2*Ml*yv[l] *Rout+M2*Ll*Vcin(t)*Rin+M2*Ll*Vcin(t ) *Rout 
-M2*Ll*yv  [5] *Rcout*Rin-*M2*Ll*yv [5] *Rcout*Rout-M2*Ll*yv [5] 
*Rcin*Rin-M2*Ll*yv [8] *Rcin*Rout ) 

/( (L2*L1*L4+L2*L1*L3-L2*M1*M1-L1*M2*M2) * (Rin+Rout ) ) ; 
dydt [5]=-(-L2*Ml*Rout*Vso (t )+L2*Ml*Rout*yv [3] *Rin+L2*Ml*Rout 

*yv[4] *Rin+L2*Ml*yv [1] *Rin+L2*Ml*yv[l] *Rout-L2*Ll*Vcin(t) 
*Rin-L2*Ll*Vcin(t)*Rout+L2*Ll*yv [5] *Rcout*Rin+L2*Ll*yv [5] 
*Rcout*Rout+L2*Ll*yv [5] *Rcin*Rin+L2*Ll*yv [5] *Rcin*Rout+Rout 
*Vso(t)*Ll*M2-Rout*yv[3] *Rin*Ll*M2-Rout*yv [4] *Rin*Ll*M2 
-yv[2] *Rin*Ll*M2-yv [2] *Rout*Ll*M2)/( (L2*L1*L4+ 
L2*L1*L3-L2*M1*M1-L1*M2*M2) * (Rin+Rout ) ) ; 
dydt  [6]  =  Po*yv[l]; 


> 
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int  main(){ 

FILE  *f ile; 

float  *vstart  =  vector(l ,N0VAR) ; 
float  result; 
int  i ; 

int  nok,  nbad; 

/♦Now  we  wish  to  merely  step  through  the  algorithm  and  determine 

a  set  of  results.  We  will  be  using  just  a  fourth  order  Runge-Kutta 
method  as  of  now.  I  may  change  this  later.  We  are  merely 
solving  the  differential  equations  for  a  JJ  pair  assumed  to  be 
in  the  configuration  of  a  squid.  We  can  then  view  the  results 
which  will  be  placed  in  a  file  for  viewing 


The  following  configurations  will  be  used. 

We  shall  drive  a  current  through  the  system  that  is  just  under  the 
maximum  critical  current  of  2.0e-4  amps.  We  will  then  apply  a 
control  current  to  the  control  circuitry  linked  to  our  JJ>s 
through  the  inductance  L.  This  should  alter  the  configuration 
of  the  system  and  create  a  HVS .  We  shall  analyze  the  input  current, 
the  currents  through  each  JJ,  and  the  voltage  across  each  JJ 
to  determine  what  state  the  system  is  in! 

Later  models  will  add  the  output  resistance  to  shunt  the  current 
along  with  an  input  current  through  a  voltage  and  a  resistor.  An 
inductance  will  be  used  later  to  model  effects  of  the  loop  and  its 
self  inductance. 

We  want  to  see  what  the  basic  operation  of  the  SQUID  should  be*/ 


xp  =  vector (l,kmax) ; 

yp  =  matrix (1 ,N0VAR, 1 ,kmax) ; 

/♦Clear  the  Matrices*/ 

f or ( i=l ; i<=N0VAR ; i++) { 
vstart [i]=0.0; 

> 

printf ("Performing  Runge  Kutta  5th  order  adaptive  integration. \n") ; 
odeint (vstart ,N0VAR ,0.0 , SIMTIM , EPS , 1 . Oe-13 , 1 . Oe- 16 , &nok , 

&nbad , derive , rkqs ) ; 

printf  ("nok=  '/.d  nbad=  */,d\n"  , nok, nbad) ; 


printf ("Simulation  over,  calculating  and  printing  results\n"); 
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/♦We  now  should  have  our  answer  and  we  just  need  to  print  out*/ 
file  =  f open("JJCONj j " , Mw") ; 
f or(i=l ; i<=kount ; i++) 

fprintf  (file  /"/,g  */,g\nM  ,xp [i]  ,yp[5]  [i]  ) ; 
f close(f ile) ; 

/♦Print  out  the  input  current  */ 
file  =  fopen("JJINj  j"  /  V) ; 
f or(i=l ; i<=kount ; i++){ 

result=(Vso(xp[i] )-(Vso(xp[i] )*Rout-Rin*Rout*(yp [3]  [i] 
+yp[4]  [i] ) )/(Rin+Rout) )/Rin; 
fprintf  (f  ile ,  '"/,g  ,/<g\n"  ,  xp  [i]  ,  result ) ; 

} 

f close(f ile) ; 

/♦Print  out  the  output  current*/ 
file  =  f open("JJOUTj j  M , Mw" ) ; 
f  or ( i= 1 ; i <=kount ; i++ ) { 

result=((Vso(xp[i] )*Rout-Rin*Rout*(yp[3] [i]+yp[4] [i] )) 

/ (Rin+Rout ) ) /Rout ; 

fprintf  (file/"/, g  */,g\n"  ,xp[i]  ,  result) ; 

} 

f close(f ile) ; 

/♦Print  out  II  */ 

file  =  f open(M JJI1 j j " , "w" ) ; 

f  or ( i=0 ; i<=kount ; i++ ) 

fprintf  (file/"/, g  '/,g\nM  ,xp[i]  ,yp[3]  [i]  ) ; 
f close(f ile) ; 

/♦Print  out  12*/ 

file  =  f open(" JJI2j j " , "w" ) ; 

for(i=0; i<=kount ; i++) 

fprintf  (file/"/, g  y.gXn"  ,xp[i]  ,yp[4]  [i]  ) ; 
f close(f ile) ; 

/♦Print  out  Itotal*/ 
file  =  f open(" JJITj j " /'w”) ; 
f or(i=0 ; i<=kount ; i++) 

fprintf  (file/"/, g  y,g\n"  ,xp[i]  ,yp[3]  [i]+yp[4]  [i]  ) ; 
f close(f ile) ; 

/♦Print  out  the  total  voltage  */ 
file  =  f open(" JJVOLj j " , MwM) ; 
f or(i=l ; i<=kount ; i++) 

fprintf  (file ,  '"/,g  y,g\n"  ,  xp  [i]  ,  Vso  (xp  [i]  ) ) ; 
f close(f ile) ; 

printf  ("Done . \n") ; 


exit (0) ; 


Chapter  4 

MVTL  Simulation  Data 


The  section  contains  the  FD-TLM  Data  Set,  the  FastHenry  extraction  data  set, 
MAPLE  V  results,  and  the  C  program  for  performing  the  conventional  circuit  simu¬ 
lation  for  the  MVTL  circuit. 


4.1  FD-TLM  Data  Set 

NMVTL 

T  JJ  MVTL  Ciruit 

♦Generated  by  FDTGRAPH,  copyright  1993,  by  Christopher  G.  Sentelle 
♦Data  in  format  to  be  used  by  FD-TLM  copyright  by  Dr.  Robert  H.  Voelker 
♦University  of  Nebraska-Lincoln 
* 

* 

♦Created  on:  Fri  Aug  11  12:40:33  1995 
♦Medium  material  used  throughout 
♦Relative  permittivity 
E  1  50  1  50  1  50  1  1  1 

♦Conductivity  throughout 
L  1  50  1  50  1  50  0  0  0 
♦Magnetic  susceptibility 
M  1  50  1  50  1  50  0  0  0 
♦Relative  permeability 
U  1  50  1  50  1  50  1  1  1 
♦Add  an  insulator  of  Si02 
E  1  50  131  50  3. 53. 53. 5 
♦Substrate  Materials 

*  Josephson  Junction  #1 

R  17  19  3  17  22  3.039e+15  0.0001  0.004  0  0.1  0.002  0.0001  5e-13  5e-ll 

♦  Josephson  Junction  #2 

R  17  19  3  38  43  3.039e+15  0.0003  0.012  0  0.3  0.002  0.0001  1.5e-12  5e-ll 
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*  Josephson  Junction  #3 

R  23  26  5  18  21  3.039e+15  0.0001  0.004  0  0.1  0.002  0.0001  5e-13  5e-ll 

*  Coupling  H  field 
H  13  27  4  23  37  1  2 

*3  dimensional  object,  Superconductor 
L  28  33  4  4  17  43  -1  -1  -1 
♦Top  Surface 

L  28  33  5  5  17  43  -1  -2  -1 

L  34  34  5  S  17  43  -2  -2  -1 

L  28  33  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 

L  28  33  4  4  44  44  -1  -1  -2 

L  34  34  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 

L  34  34  4  4  17  43  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  21  27  4  4  17  22  -1  -1  -1 

♦Top  Surface 

L  21  27  5  5  17  22  -1  -2  -1 

L  28  28  S  5  17  22  -2  -2  -1 

L  21  27  5  5  23  23  -1  -2  -2 

♦Right  Side  Surface 

L  21  27  4  4  23  23  -1  -1  -2 

L  28  28  4  4  23  23  -2  -1  -2 

♦Back  Side  Surface 

L  28  28  4  4  17  22  -2  -1  -1 

*3  dimensional  object,  Superconductor 

L  21  27  4  4  38  43  -1  -1  -1 

♦Top  Surface 

L  21  27  5  5  38  43  -1  -2  -1 

L  28  28  5  5  38  43  -2  -2  -1 

L  21  27  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 

L  21  27  4  4  44  44  -1  -1  “2 

L  28  28  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 

L  28  28  4  4  38  43  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  13  19  4  4  17  22  -1  -1  -1 

♦Top  Surface 

L  13  19  5  5  17  22  -1  -2  -1 

L  20  20  5  5  17  22  -2  -2  -1 

L  13  19  5  5  23  23  -1  -2  -2 

♦Right  Side  Surface 

L  13  19  4  4  23  23  -1  -1  -2 

L  20  20  4  4  23  23  -2  -1  -2 

♦Back  Side  Surface 

L  20  20  4  4  17  22  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  13  19  4  4  38  43  -1  -1  -1 

♦Top  Surface 
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L  13  19  5  5  38  43  -1  -2  -1 

L  20  20  5  S  38  43  -2  -2  -1 

L  13  19  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 

L  13  19  4  4  44  44  -1  -1  -2 

L  20  20  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 

L  20  20  4  4  38  43  -2  -1  -1 

♦3  dimensional  object.  Superconductor 

L  7  12  4  4  17  43  -1  -1  -1 

♦Top  Surface 

L  7  12  5  5  17  43  -1  -2  -1 

L  13  13  5  5  17  43  -2  -2  -1 

L  7  12  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 

L  7  12  4  4  44  44  -1  “1  -2 

L  13  13  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 

L  13  13  4  4  17  43  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  21  21  3  3  17  22  -1  -1  -1 

♦Top  Surface 

L  21  21  4  4  17  22  -1  -2  -1 

L  22  22  4  4  17  22  -2  -2  -1 

L  21  21  4  4  23  23  -1  -2  -2 

♦Right  Side  Surface 

L  21  21  3  3  23  23  -1  -1  -2 

L  22  22  3  3  23  23  -2  -1  -2 

♦Back  Side  Surface 

L  22  22  3  3  17  22  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  21  21  3  3  38  43  -1  -1  -1 

♦Top  Surface 

L  21  21  4  4  38  43  “1  -2  -1 

L  22  22  4  4  38  43  -2  -2  -1 

L  21  21  4  4  44  44  -1  -2-2 

♦Right  Side  Surface 

L  21  21  3  3  44  44  -1  -1  -2 

L  22  22  3  3  44  44  -2  -1  -2 

♦Back  Side  Surface 

L  22  22  3  3  38  43  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  17  21  2  2  17  22  -1  -1  -1 

♦Top  Surface 

L  17  21  3  3  17  22  -1  -2  -1 

L  22  22  3  3  17  22  -2  -2  -1 

L  17  21  3  3  23  23  -1  -2  -2 

♦Right  Side  Surface 
L  17  21  2  2  23  23  -1  -1  -2 

L  22  22  2  2  23  23  -2  -1  -2 

♦Back  Side  Surface 
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L  22  22  2  2  17  22  -2  -1  -1 
*3  dimensional  object,  Superconductor 


L 

17  21  2 

2 

38 

43  -1  -1  -1 

♦Top  Surface 

L 

17  21  3 

3 

38 

43  -1  -2  -1 

L 

22  22  3 

3 

38 

43  -2  -2  -1 

L 

17  21  3 

3 

44 

44  -1  -2  -2 

♦Right  Side  Surface 

L 

17  21  2 

2 

44 

! 

t-*- 

1 

M- 

I 

to 

L 

22  22  2 

2 

44 

44  -2  -1  -2 

♦Back  Side 

Surface 

L 

22  22  2 

2 

38 

43  -2  -1  -1 

♦Infinitely  Thin  object,  xz-plane  DAMP  RES  1  OHM 
L  25  26  4  4  23  37  -2  -2  5.0e+06 

L  27  27  4  4  23  37  -2  -2  5.0e+06 

*3  dimensional  object,  Superconductor 
L  34  39  4  4  38  43  -1  “1  -1 
♦Top  Surface 

L  34  39  5  5  38  43  -1  -2  -1 

L  40  40  5  5  38  43  -2  -2  -1 

L  34  39  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 

L  34  39  4  4  44  44  -1  -1  -2 

L  40  40  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 

L  40  40  4  4  38  43  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  40  45  4  4  3  43  -1  -1  -1 

♦Top  Surface 

L  40  45  5  5  3  43  -1  -2  -1 

L  46  46  5  5  3  43  -2  -2  -1 

L  40  45  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 
L  40  45  4  4  44  44  -1  -1  -2 

L  46  46  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 
L  46  46  4  4  3  43  -2  -1  -1 

♦Infinitely  Thin  object,  xz-plane  SRC  RES  60  OHMS 

L  40  45  4  4  2  2  -2  -2  2380.95 

L  46  46  4  4  2  2  -2  -2  2380.95 

♦3  dimensional  object,  Superconductor 

L  46  49  4  4  38  43  -1  -1  -1 

♦Top  Surface 

L  46  49  5  5  38  43  -1  -2  -1 

L  50  50  5  5  38  43  -2  -2  -1 

L  46  49  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 
L  46  49  4  4  44  44  -1  -1  -2 

L  50  50  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 
L  50  50  4  4  38  43  -2  -1  -1 
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♦Infinitely  Thin  object,  xz-plane  LOAD  RES  16  OHM 
L  50  50  4  4  38  43  8928.57  -2  -2 
L  50  50  4  4  44  44  8928.57  -2  -2 
*3  dimensional  object,  Superconductor 
L  28  33  6  6  3  47  -1  -1  -1 
♦Top  Surface 

L  28  33  7  7  3  47  -1  -2  -1 

L  34  34  7  7  3  47  -2  -2  -1 

L  28  33  7  7  48  48  -1  -2  -2 

♦Right  Side  Surface 

L  28  33  6  6  48  48  -1  -1  -2 

L  34  34  6  6  48  48  -2  -1  -2 

♦Back  Side  Surface 

L  34  34  6  6  3  47  -2  -1  -1 

*3  dimensional  object.  Superconductor 

L  2  27  6  6  46  47  -1  -1  -1 

♦Top  Surface 

L  2  27  7  7  46  47  -1  -2  -1 

L  28  28  7  7  46  47  -2  -2  -1 

L  2  27  7  7  48  48  -1  -2  -2 

♦Right  Side  Surface 

L  2  27  6  6  48  48  -1  -1  -2 

L  28  28  6  6  48  48  -2  -1  -2 

♦Back  Side  Surface 

L  28  28  6  6  46  47  -2  -1  -1 

♦3  dimensional  object.  Superconductor 

L  2  3  6  6  6  45  -1  -1  -1 

♦Top  Surface 

L  2  3  7  7  6  45  -1  -2  -1 

L  4  4  7  7  6  45  -2  -2  -1 

L  2  3  7  7  46  46  -1  -2  -2 

♦Right  Side  Surface 

L  2  3  6  6  46  46  -1  -1  -2 

L  4  4  6  6  46  46  -2  -1  -2 

♦Back  Side  Surface 

L  4  4  6  6  6  45  -2  -1  -1 

*3  dimensional  object,  Superconductor 

L  4  26  6  6  6  9  -1  “1  -1 

♦Top  Surface 

L  4  26  7  7  6  9  -1  -2  -1 

L  27  27  7  7  6  9  -2  -2  -1 

L  4  26  7  7  10  10  -1  -2  -2 

♦Right  Side  Surface 

L  4  26  6  6  10  10  -1  -1  -2 

L  27  27  6  6  10  10  -2  -1  -2 

♦Back  Side  Surface 

L  27  27  6  6  6  9  ”2  -1  -1 

*3  dimensional  object,  Superconductor 

L  24  26  6  6  10  17  -1  -1  -1 

♦Top  Surface 

L  24  26  7  7  10  17  -1  -2  -1 
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L  27  27  7  7  10  17  -2  -2  -1 

L  24  26  7  7  18  18  -1  -2  -2 

♦Right  Side  Surface 

L  24  26  6  6  18  18  -1  -1  -2 

L  27  27  6  6  18  18  -2  “1  -2 

♦Back  Side  Surface 

L  27  27  6  6  10  17  -2  -1  -1 

♦3  dimensional  object,  Superconductor 

L  22  26  6  6  18  22  -1  -1  -1 

♦Top  Surface 

L  22  26  7  7  18  22  -1  -2  -1 

L  27  27  7  7  18  22  -2  -2  -1 

L  22  26  7  7  23  23  -1  -2  -2 

♦Right  Side  Surface 
L  22  26  6  6  23  23  -1  -1  -2 

L  27  27  6  6  23  23  -2  -1  -2 

♦Back  Side  Surface 
L  27  27  6  6  18  22  -2  -1  -1 

♦Infinitely  Thin  object,  xz-plane  CNTR  SRC  RES  16  OHMS 
L  28  33  6  6  2  2  -2  -2  8928.87 

L  34  34  6  6  2  2  -2  ”2  8928.57 

♦Infinitely  Thin  object,  xz-plane  CNTR  LOAD  RES  1  OHM 

L  1  1  6  6  45  47  250000  -2  -2 

L  1  1  6  6  48  48  250000  -2  -2 

♦3  dimensional  object.  Superconductor 

L  1  6  4  4  38  43  -1  -1  -1 

♦Top  Surface 

L  1  6  5  5  38  43  -1  -2  -1 

L  7  7  5  5  38  43  -2  -2  -1 

L  1  6  5  5  44  44  -1  -2  -2 

♦Right  Side  Surface 
L  1  6  4  4  44  44  -1  -1  -2 

L  7  7  4  4  44  44  -2  -1  -2 

♦Back  Side  Surface 
L  7  7  4  4  38  43  -2  -1  -1 

♦  1  ps  =  600  iterations 

♦ 

*  Minimum  Grid  Spacing 
A  le-06 

* 

♦Simulation  time  (ps)  120 
S  111000 

♦Backup  Interval  (ps)  125 
B  120000 

♦Plot  Interval  (ps)  0.1 
P  60 
* 

♦  A  Voltage  Source 

♦  A  Pulse  waveform 

♦  Zero  Initial  Time  0(ps) 

♦  Rise  Time  5(ps) 


*  On  Time  lOO(ps) 

*  Fall  Time  5(ps) 

*  On  Voltage  -0.02 

*  Off  Voltage  0 

* 

VP  40  46  441  1Z0  3000  60000  3000  -0.02  0 

*  A  Voltage  Source 

*  A  Pulse  waveform 

*  Zero  Initial  Time  40(ps) 

*  Rise  Time  5(ps) 

*  On  Time  20(ps) 

*  Fall  Time  5(ps) 

*  On  Voltage  -0.004 

*  Off  Voltage  0 

* 

V  P  28  34  6  6  1  1  Z  24000  3000  12000  3000  -0.004  0 
*Voltage  Paths 
W  Y  1  3  42  6 
W  Y  1  5  31  8 
W  Y  1  5  20  46 
W  Y  1  3  27  19 
W  Y  1  3  27  41 
W  Y  1  3  14  20 
W  Y  1  3  13  40 
W  Y  1  3  36  39 
W  Y  1  3  46  41 
* 

^Current  Loops 
J  Z  39  46  3  5  6 
J  Z  27  34  5  7  8 
J  X  5  7  45  48  20 

J  X  3  5  16  23  27 

J  X  3  5  37  44  27 

J  X  3  5  16  23  14 

J  X  3  5  37  44  13 

J  X  3  5  37  44  36 

J  X  3  5  37  44  47 

*  Variable  Mesh  Array 

*  Variable  Mesh  in  the  X  direction 
G  X  1  50  1 

*  Variable  Mesh  in  the  Y  direction 
G  Y  1  50  1 

*  Variable  Mesh  in  the  Z  direction 
G  Z  1  50  1 
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4.2  FastHenry  Data  Set 


This  is  the  data  set  provided  to  the  FastHenry  program  to  determine  the  values  of 
the  two  parasitic  inductances  in  the  MVTL  circuit. 

♦Determination  of  Inductance  for  Squid  Loop 

♦  June  2,  1995  ,  Christopher  Sentelle 

♦  Uses  fasthenry 

♦  This  is  the  inductances  used  for  fmvtl4.dat 

♦  where  the  inductance  values  had  been  updated. 

.Units  um 

♦  Use  a  high  value  for  conductivity,  Superconductor 
.Default  nhinc  =  1  nwinc  =  5  sigma=1.0e20  z=3  w=6  h=l 

♦  Ground  Plane 

gl  xl=0  y 1=0  zl=0 
+  x2=50  y2=0  z2=0 

+  x3=50  y3=50  z3=0 

+  segl=20  seg2=20 
+  thick=l 

♦  Setup  for  the  system 

♦  Squid  Loop 
N1  x=19  y=36 
N2  x=19  y=16 
N3  x=22  y=19 
N4  x=3T  y=19 
N5  x=40  y=16 
N6  x=40  y=36 
N7  x=40  y=37 
N8  x=40  y=43 
N9  x=37  y=40 
N10  x=22  y=40 
Nil  x=19  y=43 
N12  x=19  y=37 
El  Ml  N2 

E2  N3  N4 
E3  N5  N6 
E4  N7  N8 
E5  N9  N10 
E6  Nil  N12 
.equiv  N2  N3 
.equiv  N4  N5 
.equiv  N8  N9 
.equiv  N10  Nil 
♦Control  line 
N13  x=19  y=2  z=5 
N14  x=19  y=47  z=5 
N15  x=22  y=46  z=5 
N16  x=46  y=46  z=5 
N17  x=47  y=47  z=5 
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N18  x=47  y=5  z=5 

N19  x=46  y=7  z=5 

N20  x=26  y=7  z=5 

N21  x=24 . 5  y=5  z=5 

N22  x=24 . 5  y=17  z=5 

E7  N13  N14 

E8  N15  N16  w=2 

E9  N17  N18  w=2 

E10  N19  N20  w=4 

Ell  N21  N22  w=3 

*  Used  to  determine  Mf 

E12  N6  N7 

.equiv  N14  N1S 

. equiv  N16  N17 

.equiv  N18  N19 

.equiv  N20  N21 

♦Define  External  Ports 

.external  N1  N12 

.external  N13  N22 

.freq  fmin=le9  fmax=le9  ndec=l 

.  end 
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>  eqnl  :=Icl*sin( phil ( t ) ) +Cjl*dif f (VI ( t ) . t ) +G1 (VI ( t ) )  +  (Vl( t ) -V2( t ) ) 

>  /Rd-Il ( t )+Ic3*sin(phi3 ( t ) )+G3 ( VI ( t )-V4 ( t ) )+Cj 3*di f f ( VI ( t ) -V4 ( t ) , 

>  t)=0; 

eqnl  =  Id  sin(phil(/))  +  Cjl  Vl(oj  +  Gl(  V1(0)  +  ~  *U0 

+  /c5sin(phi3(O)  +  G3(Vl(O-V4(/))  +  Q5[[|vi(/)]-^|v4(/)jj  =  0 

>  eqn2 : =Ic2*sin(phil ( t )-2*PI* (Mf *13 ( t )+L2*I2 ( t ) -L1*I1 ( t ) )/Flo)+Cj2 

>  *dif f (V2 ( t ) , t )+G2 (V2( t ) )+(V2( t )-Vl ( t ) )/Rd-I2 ( t )=0 : 


eqn2  :=  lc2  soil  phi  1  ( ()  -  2 - - - — - I  +  Cj2 


+  G2(  V2  (0)  +  -  12(f)  =  0 

>  eqn3 :=Il(t)+I2(t)+(V3(t)-Vin(t) )/Rin+I4 ( t )=0 : 

eqn3  :=  11(0  + 12(0  +  V3(/)-  V'n(0  +  I4(,}  =  0 

>  eqn4:=V4( t )/Rcout-I3 ( t )-Ic3*sin (phi3 ( t ) ) -G3 (VI ( t ) -V4 ( t ) ) -Cj 3*dif 

>  f(Vl(t)-V4(t),t)=0: 


eqn4  := 

V4(0 


-^-I3(O-/cisin(phi3(/))-G3(Vl(/)-V4(O)-Q'3^Vl(/)J-^-V4(/)JJ  =  0 

>  eqn5 : =-I4 ( t ) +V5 ( t )/Rout=0 ;  ' 

eqn5  :=  -I4(  0  +  =  0 

>  eqn6 :=I3(t)=(Vcin(t)-(L3*diff(I3(t) , t ) -Ml*dif f ( II ( t ) , t )+M2*dif f ( 

>  12 ( t ) , t ) +V4 ( t ) ) )/Rcin : 


eqn6  :=  13(0  =  " 


Vcin(/)-i3  ^-13(0  +A//  ^I1(/)  -M2  12(0  -V4(0 


>  eqn7 : =V3 ( t )-Vl ( t ) =Ll*dif f ( II ( t ) , t ) -Ml’dif f ( 13 ( t ) , t ) ; 

eqn7  :=  V3(/)  -  Vl(/)  =  LI  Il(o)  -  Ml  I3(o) 

>  eqn8:=V3(t)-V2(t)=L2*diff (I2(t) . t ) +M2*dif f ( 13 ( t ) .  t); 


eqn8  :=  V3(0  -  V2(0  =  L2  hr  12(0  +  M2  hr  13(0 


>  eqn9:=diff (phil(t) , t)=Po*Vl(t) : 


eqn9  :=^phil(/)  =  Po  VI (/) 

>  eqnlO : =di f f ( phi2 (t ) . t ) =Po*V2 ( t ) : 

eqnl 0  phi2(0  =  Po  V2(  0 

ct 

>  eqnll : =dif f (phi3 ( t } . t )=Po* { VI ( t )-V4 ( t ) ) ; 


a 

eqnll  :=  — phi3(/)  =  Po(Vl(/)-  V4(f)) 


>  eqnl2:=V3(t)-V5(t)=L4*diff (I4(t) ,t)  ; 

eqnll  :=  V3(f)  -  V5(f)  =  L4  14(f)] 
7  auxl : =solve ( eqn3 ~ V3 ( t ) ) ; 


auxl  :=  - 

[ll(f)  +  I2(f)-^i  +  I4(f)^ 

|  Rin 

>  aux2 : =solve ( eqn5 , V5 ( t ) ) ; 

aux2  ~  14 (t)  Rout 

>  eqn7a : =subs( V3 ( t )=auxl , eqn7 ) ; 

eqn7a  :=-^Il(/)  + 12(0  14(')]  Pin  -  Vl(f)  =  LI  ^  11(f)] -  Ml  ^  13(f)] 

>  eqn8a:=subs(V3(t)=auxl.eqn8);  “ 

eqn8a  :=  -^ll(f)  +  12(f)  -  +  14(f)]  Rin  -  V2(f)  =  L2  ^  12(f)]  +  M2  ^  13(f)] 

>  eqnl2a : =subs ( {V5 ( t ) =aux2 . V3 ( t ) =auxl } . eqnl2 ) ; 

eqnlla  :=  |ll(f  )  +  12(f)  -  — ^  +  14(f)]  Rin  -  14(f)  Rout  =  L4  ^  14(f)] 

>  ans : =solve( {eqnl , eqn2 , eqn4  ,  eqn6  .  eqn7a . eqn8a. eqn9 , eqnll , eqnl2a}  { 

>  diff(Vl(t).t),diff(V2(t).t}.diff(V4(t).t),diff(Il(t).t),diff(I2( 

>  tj.t) ,diff(I3(t).t) .diff(I4(t).t) .di£f(phil(t).t).diff (phi3(t).t 

am  :=  ff  phil(  =  Po  V1('X  f  Phi3(0  =  Po  (VI  (f )  -  V4(f )),  |  V4(f)  =  -  ( 

Q3  Rcoutlcl  sin(phil(f))/W+  Cj3  Rcout  Gl(  Vl(f ))  Rd-  Cj3  Rcout  V2 (/) 

+  Cj3  Rcout  Vl(f)-QJ  Rcout  U(t)Rd+Cj3  RdV4(t)  -  Cj3  /W  13(f)  Rcout 
+  Cjl  Rd~V4(t )  -  Cjl  Rdl3(t)  Rcout- Cjl  Rd Ic3  sin( phi3(f ))  Rcout 

-  Cjl  &/G3(Vl(f)- V4(t))Rcout)/(Cjl  Rd  Rcout  Cj3),j\\(t)  =  -{n  Ml  U(t)Rcin 

-  L2  Ml  Vcin(f)  +  L2  Ml  V4(f)  +  L2  L3  11(f)  Rin  +  L2 L3 12(f)  Rin  -  L2  L3  Vin(f) 

+  L2  L3  14(f)  Rin  +  L2  L3  Vl(f )  -  M2 11(f)  Rin  Ml  -  M2  \2(t)Rin  Ml  +  M2  Vin(f)  Ml 

-  M2 14(f)  Rin  Ml  -  M2  V2(f)  Ml - M2 2  II  (l)Rin-M22 12(f)  Rin  +  M22  Vin(f) 

-  M22  14(f)  Rin  -  M2 2  Vl(f  ))/(-M22  LI  +  L2  L3  LI  -  L2  Ml2),  |  I2(f )  =  -  ( 

-13 (t)Rcin  M2  LI  +  Ml 2  Vin(f)  -  Ml2  V2(f )  +  M2  Vin(f)  Ml  -  Ml  M2  Vl(f) 

+  L3  LI  11(f)  Rin  +  L3  LI  12(f)  Rin  -  L3  LI  Vin (t)+L3  LI  14(f)  Rin  +  L3  LI  V2(f) 

-  Ml2  11(f)  Rin  -  Ml 2  12(f)  Rin  -  Ml 2  14(f)  Rin  -M2  11(f)  Rin  Ml  -  M2  12(f)  Rin  Ml 

-  M2 14(f)  Rin  MI  +  Vcin(f)  M2  LI  -  V4(f)  M2  Ll)/{-M22  LI  +  L2  L3  LI  -  L2  Ml2), 

j'l3(t)  =  -(-Vir\(t)L2  Ml  +  V](t)L2Ml  +L1  L2l3(t)  Rein-  LI  Z,2Vcin(f) 


> 


+  LI L2  V4 (t)-LI  M2  II (t)Rin-Ll  M2 12(c)  Rin  +  LI  M2  Vin(/)  -  LI  M2 14(c)  Rin 
—  LI  M2  V2(c)  +  11(c)  Rin  L2  Ml  +  12(C)  Rin  L2  Ml  +  14(c)  Rin  L2  Ml ) /( 

-M2  LI  +  L2  L3  LI  -  L2 Ml 2),  —  Vl(c)  =  -  ( Rcout  Id  sin(phil(c))  Rd 


+  Rcout  Gl(  Vl(/))  Rd- Rcout  V2(C)  +  Rcout  VI (C)  -  Rcout  11(c)  Rd+RdV4(t) 

-  Rd  13(c)  Rcout )/ (  Cjl  Rd  Rcout),  ^  V2(c)  =  ^ 

.J^,Jvmt)Flo-2YlMfYi(t)-2UL2Y2(t)  +  2ULimt)^^. 

“  Sl\  Jfo - — J  Rd~  G2(  V2(c))  /W 

-  V2(c)  +  Vl(c)  +  12(c)  «rfj/(Q2  Rd), 

J-I4(C)  =  I1(/)  —  + 12(/)  Rin  ~  Vin(0  + 14( C )  Rin  +  14(C)  Rom | 


ans[l] ; 


|phil(c)=p0vi(c) 

>  ans[2 ]  ;  ~  ~  — — : — - - — 

|phi3(c)  =  /’0(Vl(c)-V4(c)) 

>  an s  [  3  ] ;  ~~  — - - - = - 

d 

-  V4(c)  =  -(CJ3  Rcout  Id  sin( phi  1(c))  Rd+  Cj3  Rcout  Gl(  Vl(c))  Rd-  Cj3  Rcout  V2(c) 

+  CJ3  Rcout  VI (c)  -  Cj3  Rcout  II {t)Rd+  Cj3  RdVA{t)  -  Cj3  RdY3{t)  Rcout 
+  Cjl  Rd\4(t)  -  Cjl  Rd  13(c)  Rcout  -  Cjl  Rd Ic3  sin(phi3(c))  Rcout 
_ Z£E  RdG3(  Vl(t )  -  V4( C ) )  Rcout )/(Cjl  Rd  Rcout  Cj3 ) 

>  ans[ 4]  ;  '  “  - - - - - - - 

Jt  11(c)  =  -  ( L2  Ml  13(C)  Rein  - L2 Ml  Vein (c)  +  L2  Ml  V4(c)  +  L2  L3 11(c)  Rin 
+  L2  L3  12(C)  Rin  -  L2  L3  Vin(C)  +  L2  L3 14(c)  Rin  +  L2  L3  Vl(c)  -  M2 11(C)  Rin  Ml 

-  M2 12(c)  Rin  Ml  +  M2  Vin(c)  Ml  -M2 14(c)  Rin  Ml  -M2  V2(c)  Ml  -M21 11(c)  Rin 

-  M2 2 12(C)  Rin  +  M22  Vin(c)  -  M22 14(c)  Rin  -  M22  Vl(c)) /( 

-M22  LI  +  L2  L3  LI  -  L2  Ml 2) 

>  ans[5 ]  ;  "  - - - - — - - - 

^  12(c)  =  -  (-13(c)  Rein  M2  LI  +  Ml 2  Vin(c)  -  Ml2  V2(c)  +  M2  Vin(c)  Ml  -  Ml  M2  VI (c) 
+  L3  LI  11(c)  Rtn  +  L3  LI  \2{t) Rin  - L3  LI  Vin (t)  +  L3Ll  14 (t)Rin  +  L3 LI  V2(C) 
-Ml2 11(c)  Rin -Ml2  l2(t)Rin  -Ml2 14(c)  Rin-M2\\(t)  Rin  MI-M2Y2{t)  Rin  Ml 

-  M2 14(c)  Rin  Ml  +  Vcin(c)  M2  LI-  V4(c)  M2  Ll)/(-M22  LI  +  L2  L3  LI  -  L2  Ml2) 

y~ans[sY;  '  ~ — - - - — - - - 
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— 13(/)  =  -  (-Vin(/)  L2  Ml  +  VI (f)  L2  Ml  +  LI  L2 13(f)  Rein  -  LI  L2  Vein (t)  +  Ll  L2  V4(t) 

-  LI  M2 11(f)  Rin  -LI  M2 12(f)  Rin  +  LI  M2  Vin(f )  -LI  M2 14(  t )  Rin  -LI  M2  V2(f ) 

+  11(0  Rm  L2  Ml  +  12(f)  Rin  L2  Ml  +  14(f)  Rin  L2  Ml )/{-M22  LI  +L2  L3  LI  -L2  Ml2) 

>  ans[7] ;  "  ;  “ 

—  Vl(f)  =  -  ( Rcout  lei  sin( phi  1(f))  Rd  +  Rcout  Gl(  Vl(f))  Rd-Rcout  V2(f)  +  Rcout  VI (f) 

-  Rcout  l\(t)Rd+Rd  V4(f )  -  Rd  13(f)  Rcout  )f(Cjl  Rd  Rcout) 

>  ans[8]:  ““  ~  —■ — — —  —  — 

|  V2(f)  =  ^-Ic2  :in^phil(OF/t>~2nA//I3(°~2nZ2I2(<)  +  2nZ/I1(f)j  Rd 


-  G2(  V2(f))  Rd-  V2(f)  +  Vl(f)  +  12(f)  Rd ]/( Cj2  Rd) 


>  ans[9j; 


-Mt)=. 


Il(  f)  Rin  + 12(  f )  Rin  -  Vin(f )  +  14(f)  Rin  +  14(f)  Rout 
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4.4  C  Code 


This  is  a  portion  of  the  C  Code  used  to  solve  the  set  of  first  order  differential  equations 

provided  by  MAPLE  using  the  Runge-Kutta  fifth-order  method.  Conventional  circuit 

simulation  results  are  then  obtained. 

#include  "nrutil.c" 

#include  "odeint.c" 

#include  "rkqs.c" 

#include  “rkck.c" 


#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 
#def ine 


NOVAR  6 
L3  6 . Oe-12 
L4  6. Oe-12 
Ml  2. Oe-12 
M2  2. Oe-12 
Mf  2 . 068e-12 
LI  8 . 39e-12 
L2  8 . 39e-12 
Ic  1.0e-6 
G1  4 . Oe-3 
II  0.0 
G2  0.1 
Vs  2. Oe-3 
Vt  1.0e-4 
Io  1.0e-4 
Po  3 . 039el5 
Cj  5 . Oe-13 
Rein  17.78 
Rcout  1 . 19 
Rin  71.99 
Rout  20.0 
FLo  2 . 068e-15 

TS  1.0e-15  /*0ne  femtosecond  shall  be  used*/ 

PLT  0.1e-12 

PI  3.141592763 

SIMTIM  180. Oe-12 

EPS  1.0e-2 


int  kmax=1300 ,kount ; 

float  *xp>**yp,dxsav=l. Oe-13; 


float  G(float  v) 


float  ans; 


ans  =  Gl*v  +  (II  +  G2*fabs(v))*((l/(l+exp((Vs-v)/Vt))) 

-(1/ (l+exp((Vs+v)/Vt) ) ) ) ; 

return  ans ; 

> 


float  Vso(float  t) 

{ 

float  starttime  =  0.0e-12; 
float  risetime  =  B.Oe-12; 
float  ontime  =  lOO.Oe-12; 
float  falltime  -  S.0e-12; 
float  hicur  =  0.010; 

if  (t  <=starttime) 
return  0.0; 

if  ( (t>starttime)  kk  (t<=starttime+risetime) ) 
return  (t-st arttime) *hicur/ risetime ; 
if  ( (t>starttime+risetime)&&(t<=starttime+risetime+ontime) ) 
return  hicur; 

if ( ( t >startt ime+r iset ime+ont ime ) &&(t <=st artt ime+r iset ime+ont ime 
+f alltime)) 

return  (hicur  -  (t-(startt ime+risetime+ontime) )*hicur 
/falltime); 

else 

return  0.0; 

} 

float  Vcin(float  t) 

{ 

float  starttime  =  30.0e-12; 
float  risetime  =  5.0e-12; 
float  ontime  =  20.0e-12; 
float  falltime  =  5.0e-12; 
float  hicur  =  0.003; 


if  (t  <=starttime) 
return  0.0; 

if  ( (t>starttime)  kk  (t<=starttime+risetime) ) 
return  (t-starttime)*hicur/risetime; 
if  ( (t>starttime+risetime)&&(t<=starttime+riset ime+ont ime) ) 
return  hicur; 

if ( ( t> st artt ime+r iset ime+ont ime )&&(t<=startt ime+r iset ime+ont ime 
+f alltime)) 

return  (hicur  -  (t-(starttime+risetime+ontime) ) *hicur 
/falltime) ; 


else 


} 


return  0.0; 
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void  derive  (float  t,  float  yv[]  ,  float  dydt[]) 

{ 

int  i ; 

/*In  this  new  system,  the  following  variables  apply: 
y  [1]=V1 
y  [2]  =V2 
y [3] =V4 
y  [4]  =1 1 
y[5]=I2 
y[6]=I3 
y [7] =phil 
y [8]=phi3 

*/ 

float  G(float  v); 

float  Vso(float  t); 

float  Vcin(float  t); 

dydt [l]  =  (Rd*yv[6] *Rcout-Rd*G(yv[3]-yv[l] )*Rcout-2 . *Rd*yv [6] 

*s in (yv [8] ) *Rcout-Rd*yv [3] +Rcout *y v [4] *Rd+Rcout*y v  [2] 
-Rcout*yv  [1] -Rcout*G (yv [1] ) *Rd-Rcout*yv [4] *sin(yv [7] ) 
*Rd-Rcout*G(yv[l]-yv[3] ) *Rd)/(Cj l*Rcout*Rd) ; 

dydt [2] =- (-y v  [5] *Rd+y v  [2] -yv [l] +G (yv [2] ) *Rd+Ic2*s in( (y v [7] *Flo-2 
*PI* (Mf *yv  [6] -L2*yv [S] +Ll*yv [4] ) ) /Flo ) *Rd ) / (Cj 2*Rd) ; 

dydt [3]=(Cj3*Rd*I3(t)*Rcout-Cj3*Rd*G(yv [3] -yv[l] ) *Rcout 

-2*Cj3*Rd*Ic3*sin(yv [8] ) *Rcout-Cj3*Rd*yv [3]+Cj  3*Rcout 
*yv[4] *Rd+Cj3*Rcout*yv [2] -Cj3*Rcout*yv [l] -Cj3*Rcout 
*G(yv[l] )*Rd-Cj3*Rcout*Icl*sin(yv [7] )*Rd 
-C j  3*Rcout*G (yv [1] -yv [3] ) *Rd+Cj l*Rd*yv [6]  +Rcout-C j l*Rd 
*G (yv [3] -y v  [1] ) *Rcout-C j l*Rd*Ic3*s in(y v [8] ) *Rcout-C j l*Rd 
*yv  [3] ) / (C j l*Rcout *Rd*C j  3) ; 

dydt [4] =- (Ml*L2*yv  [6] *Rcin-Ml*L2*Vcin(t )+Ml*L2*yv [3] +M2*Ml*yv  [1] 


> 


int  main(){ 

FILE  *f ile; 

float  *vstart  =  vector (1 ,N0VAR) ; 
float  result; 
int  i ; 

int  nok,  nbad; 

/*  This  next  simulation  will  simulate  the  MVTL  circuit.  Simulation 
will  be  performed  with  this  system  followed  by  simulation  with  the 
FDTLM  method. 

The  following  configurations  will  be  used. 

We  shall  drive  a  current  through  the  system  that  is  just  under  the 
maximum  critical  current  of  2.0e-4  amps.  We  will  then  apply  a 


control  current  to  the  control  circuitry  linked  to  our  JJ’s 
through  the  inductance  L.  This  should  alter  the  configuration  of 
the  system  and  create  a  HVS.  We  shall  analyze  the  input  current, 
the  currents  through  each  JJ,  and  the  voltage  across  each  JJ  to 
determine  what  state  the  system  is  in! 

Later  models  will  add  the  output  resistance  to  shunt  the 
current  along  with  an  input  current  through  a  voltage  and  a 
resistor.  An  inductance  will  be  used  later  to  model  effects 
of  the  loop  and  its  self  inductance. 

We  want  to  see  what  the  basic  operation  of  the  SQUID  should  be*/ 


xp  =  vector(l ,kmax) ; 

yp  =  matrix (1 ,N0VAR, 1 ,kmax) ; 

/♦Clear  the  Matrices*/ 

f or (i=l ; i<=N0VAR; i++){ 
vstart [i]=0. 0 ; 

> 

printf ("Performing  Runge  Kutta  Sth  order  adaptive  integration. \n") ; 
odeint (vstart , NOVAR ,0.0, SIMTIM , EPS , 1 . Oe- 13 , 1 . Oe- 16 , &nok , 

&nbad , derive , rkqs ) ; 

printf ("nok=  */,d  nbad=  '/.dXn"  ,nok,nbad) ; 


printf ("Simulation  over,  calculating  and  printing  results\n"); 
/*We  now  should  have  our  answer  and  we  just  need  to  print  out*/ 
file  =  f open("JJC0N" , "w") ; 
f or(i=l ; i<=kount ; i++) 

f  printf  (file ,  M,/*g  •/#g\nH  ,xp[i]  ,yp[5]  [i]  ) ; 
f close(f ile) ; 

/♦Print  out  the  input  current  */ 
file  =  f open ("J JIN" , "w") ; 
f or(i=l ; i<=kount ; i++){ 

result=(Vso(xp[i] )-(Vso(xp[i] )*Rout-Rin*Rout*(yp [3] [i] 
+yp[4]  [i] ) )/ (Rin+Rout) )/Rin; 
f  printf  (file ,  n,/,g  y«g\n"  ,xp[i]  , result) ; 

> 

f close(f ile) ; 

/♦Print  out  the  output  current*/ 
file  =  f open(" JJOUT" , "w" ) ; 
f  or(i=l ;  i<=kount ;  i++H 

result=( (Vso(xp[i] )*Rout-Rin*Rout*(yp[3] [i]+yp[4]  [i] ) ) 

/ (Rin+Rout ) ) /Rout ; 
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fprintf  (file , "'/,g  '/.g\n",xp[i]  .result)  ; 

> 

f close(f ile) ; 

/♦Print  out  II  */ 

file  =  f open("JJIi" , "w") ; 

f or(i=0; i<=kount ; i++) 

fprintf  (file,  n,/,g  '/,g\n"  ,xp[i]  ,yp[3]  [i] ) ; 
f close(f ile) ; 

/♦Print  out  12+/ 

file  =  fopen("JJI2","w") ; 

f or(i=0; i<=kount ; i++) 

fprintf  (file,  "*/,g  '/,g\n"  ,xp[i]  ,yp[4]  [i]) ; 
f close(f ile) ; 

/♦Print  out  Itotal*/ 
file  =  f open("JJIT" , "w" ) ; 
f or(i=0; i<=kount ; i++) 

fprintf  (file,"'/,g  '/,g\n" ,xp [i]  ,yp[3]  [i]+yp[4]  [i]  ) ; 
f close(f ile) ; 

/♦Print  out  the  total  voltage  */ 
file  =  f open("JJV0L" , "w") ; 
f or(i=i ; i<=kount ; i++) 

fprintf  (file,"'/,g  '/,g\n",xp[i]  ,Vso(xp[i]  ) ) ; 
f close(f ile) ; 

pr  intf  ( "Done . \n"  )  ; 

exit (0) ; 
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